yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #04154
[Branch ~yade-dev/yade/trunk] Rev 2184: 1. Remove some ToAxisAngle and FromAxisAngle as per https://www.yade-dem.org/wiki/Wm3%E2%86%92Eigen
------------------------------------------------------------
revno: 2184
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Mon 2010-04-26 15:58:23 +0200
message:
1. Remove some ToAxisAngle and FromAxisAngle as per https://www.yade-dem.org/wiki/Wm3%E2%86%92Eigen
2. Fix NewtonIntegrator using [i] on quaternions.
2. Add docs for VTKRecorder
modified:
core/State.hpp
gui/qt3/GLViewer.cpp
lib/base/yadeWm3Extra.cpp
lib/base/yadeWm3Extra_dont_include_directly.hpp
lib/serialization/KnownFundamentalsHandler.tpp
pkg/common/Engine/GlobalEngine/SpheresFactory.cpp
pkg/common/Engine/PartialEngine/JumpChangeSe3.cpp
pkg/common/Engine/PartialEngine/RotationEngine.cpp
pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp
pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp
pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp
pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp
pkg/dem/Engine/GlobalEngine/VTKRecorder.hpp
pkg/dem/Engine/PartialEngine/Disp2DPropLoadEngine.cpp
pkg/dem/Engine/PartialEngine/KinemCNDEngine.cpp
pkg/dem/Engine/PartialEngine/KinemCNLEngine.cpp
pkg/dem/Engine/PartialEngine/KinemCNSEngine.cpp
pkg/dem/Engine/PartialEngine/KinemCTDEngine.cpp
pkg/dem/PreProcessor/CohesiveTriaxialTest.cpp
pkg/dem/PreProcessor/DirectShearCis.cpp
pkg/dem/PreProcessor/STLImporterTest.cpp
pkg/dem/PreProcessor/SimpleShear.cpp
pkg/dem/PreProcessor/TriaxialTest.cpp
pkg/dem/PreProcessor/TriaxialTestWater.cpp
py/_eudoxos.cpp
py/yadeWrapper/customConverters.cpp
py/yadeWrapper/yadeWrapper.cpp
scripts/test-spiral.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/State.hpp'
--- core/State.hpp 2010-04-25 13:18:11 +0000
+++ core/State.hpp 2010-04-26 13:58:23 +0000
@@ -46,7 +46,7 @@
//! Return displacement (current-reference position)
Vector3r displ() const {return pos-refPos;}
//! Return rotation (current-reference orientation, as Vector3r)
- Vector3r rot() const { Quaternionr relRot=refOri.conjugate()*ori; Vector3r axis; Real angle; relRot.ToAxisAngle(axis,angle); return axis*angle; }
+ Vector3r rot() const { Quaternionr relRot=refOri.conjugate()*ori; AngleAxisr aa(angleAxisFromQuat(relRot)); return aa.axis()*aa.angle(); }
// python access functions: pos and ori are references to inside Se3r and cannot be pointed to directly
Vector3r pos_get() const {return pos;}
=== modified file 'gui/qt3/GLViewer.cpp'
--- gui/qt3/GLViewer.cpp 2010-04-25 12:37:47 +0000
+++ gui/qt3/GLViewer.cpp 2010-04-26 13:58:23 +0000
@@ -608,11 +608,10 @@
for(int planeId=0; planeId<renderer->numClipPlanes; planeId++){
if(!renderer->clipPlaneActive[planeId] && planeId!=manipulatedClipPlane) continue;
glPushMatrix();
- Real angle; Vector3r axis;
const Se3r& se3=renderer->clipPlaneSe3[planeId];
- se3.orientation.ToAxisAngle(axis,angle);
+ AngleAxisr aa(angleAxisFromQuat(se3.orientation));
glTranslatef(se3.position[0],se3.position[1],se3.position[2]);
- glRotated(angle*Mathr::RAD_TO_DEG,axis[0],axis[1],axis[2]);
+ glRotated(aa.angle()*Mathr::RAD_TO_DEG,aa.axis()[0],aa.axis()[1],aa.axis()[2]);
Real cff=1;
if(!renderer->clipPlaneActive[planeId]) cff=.4;
glColor3f(max((Real)0.,cff*cos(planeId)),max((Real)0.,cff*sin(planeId)),planeId==manipulatedClipPlane); // variable colors
=== modified file 'lib/base/yadeWm3Extra.cpp'
--- lib/base/yadeWm3Extra.cpp 2010-04-25 12:37:47 +0000
+++ lib/base/yadeWm3Extra.cpp 2010-04-26 13:58:23 +0000
@@ -28,7 +28,7 @@
__attribute__((deprecated)) Vector3d operator*(const Vector3d& v, const float s){return Vector3d(s*v[0],s*v[1],s*v[2]);}*/
std::ostream & operator<< (std::ostream &os, const Vector3r &v){ return os << v[0] << " " << v[1] << " " << v[2];}
-std::ostream & operator<< (std::ostream &os, const Quaternionr &q){ Vector3r axis; Real angle; q.ToAxisAngle(axis,angle); return os<<axis<<" "<<angle;}
+std::ostream & operator<< (std::ostream &os, const Quaternionr &q){ AngleAxisr aa(angleAxisFromQuat(q)); return os<<aa.axis()<<" "<<aa.angle();}
// c * M
=== modified file 'lib/base/yadeWm3Extra_dont_include_directly.hpp'
--- lib/base/yadeWm3Extra_dont_include_directly.hpp 2010-04-25 11:58:30 +0000
+++ lib/base/yadeWm3Extra_dont_include_directly.hpp 2010-04-26 13:58:23 +0000
@@ -6,6 +6,10 @@
#include<Wm3Quaternion.h>
using namespace Wm3;
+// this will be replaced by AngleAxisr(Quaternionr) ctor once we use Eigen everywhere
+template<typename RealType>
+AngleAxis<RealType> angleAxisFromQuat(const Quaternion<RealType>& q){ AngleAxis<RealType> aa; q.ToAxisAngle(aa.axis(),aa.angle()); return aa; }
+
#include"yadeWm3_dont_include_directly.hpp"
#include<limits>
=== modified file 'lib/serialization/KnownFundamentalsHandler.tpp'
--- lib/serialization/KnownFundamentalsHandler.tpp 2010-04-25 12:37:47 +0000
+++ lib/serialization/KnownFundamentalsHandler.tpp 2010-04-26 13:58:23 +0000
@@ -329,7 +329,7 @@
axis[2] = lexical_cast_maybeNanInf<RealType>(tokens[2]);
angle = lexical_cast_maybeNanInf<RealType>(tokens[3]);
}
- tmp->FromAxisAngle(axis,angle);
+ *tmp=Quaternionr(AngleAxisr(angle,axis));
}
else if (a.type()==typeid(const vector<unsigned char>*)) // from binary stream to Type
{
@@ -354,19 +354,17 @@
string * tmpStr = any_cast<string*>(a);
Quaternion<RealType> * tmp = any_cast<Quaternion<RealType>*>(ac.getAddress());
- RealType angle;
- Vector3<RealType> axis;
- tmp->ToAxisAngle(axis,angle);
- axis.normalize();
+ AngleAxis<RealType> aa(angleAxisFromQuat(*tmp));
+ aa.axis().normalize();
*tmpStr = IOFormatManager::getCustomFundamentalOpeningBracket() +
- lexical_cast<string>(axis[0]) +
- IOFormatManager::getCustomFundamentalSeparator() +
- lexical_cast<string>(axis[1]) +
- IOFormatManager::getCustomFundamentalSeparator() +
- lexical_cast<string>(axis[2]) +
- IOFormatManager::getCustomFundamentalSeparator() +
- lexical_cast<string>(angle) +
+ lexical_cast<string>(aa.axis()[0]) +
+ IOFormatManager::getCustomFundamentalSeparator() +
+ lexical_cast<string>(aa.axis()[1]) +
+ IOFormatManager::getCustomFundamentalSeparator() +
+ lexical_cast<string>(aa.axis()[2]) +
+ IOFormatManager::getCustomFundamentalSeparator() +
+ lexical_cast<string>(aa.angle()) +
IOFormatManager::getCustomFundamentalClosingBracket();
}
else if (a.type()==typeid(vector<unsigned char>*)) // from Vector2<RealType> to binary stream
@@ -429,7 +427,7 @@
axis[2] = lexical_cast_maybeNanInf<RealType>(tokens[5]);
angle = lexical_cast_maybeNanInf<RealType>(tokens[6]);
}
- orientation.FromAxisAngle(axis,angle);
+ orientation=Quaternionr(AngleAxisr(angle,axis));
*tmp = Se3<RealType>(position,orientation);
}
else if (a.type()==typeid(const vector<unsigned char>*)) // from binary stream to Type
@@ -458,12 +456,10 @@
string * tmpStr = any_cast<string*>(a);
Se3<RealType> * tmp = any_cast<Se3<RealType>*>(ac.getAddress());
- RealType angle;
- Vector3<RealType> axis;
Vector3<RealType> position;
- tmp->orientation.ToAxisAngle(axis,angle);
- axis.normalize();
+ AngleAxis<RealType> aa(angleAxisFromQuat(tmp->orientation));
+ aa.axis().normalize();
position = tmp->position;
//*tmpStr = IOFormatManager::getCustomFundamentalOpeningBracket();
@@ -476,13 +472,13 @@
//IOFormatManager::getCustomFundamentalClosingBracket() +
IOFormatManager::getCustomFundamentalSeparator() +
//IOFormatManager::getCustomFundamentalOpeningBracket() +
- lexical_cast<string>(axis[0]) +
- IOFormatManager::getCustomFundamentalSeparator() +
- lexical_cast<string>(axis[1]) +
- IOFormatManager::getCustomFundamentalSeparator() +
- lexical_cast<string>(axis[2]) +
- IOFormatManager::getCustomFundamentalSeparator() +
- lexical_cast<string>(angle) +
+ lexical_cast<string>(aa.axis()[0]) +
+ IOFormatManager::getCustomFundamentalSeparator() +
+ lexical_cast<string>(aa.axis()[1]) +
+ IOFormatManager::getCustomFundamentalSeparator() +
+ lexical_cast<string>(aa.axis()[2]) +
+ IOFormatManager::getCustomFundamentalSeparator() +
+ lexical_cast<string>(aa.angle()) +
IOFormatManager::getCustomFundamentalClosingBracket();
// IOFormatManager::getCustomFundamentalClosingBracket();
}
=== modified file 'pkg/common/Engine/GlobalEngine/SpheresFactory.cpp'
--- pkg/common/Engine/GlobalEngine/SpheresFactory.cpp 2010-03-20 12:40:44 +0000
+++ pkg/common/Engine/GlobalEngine/SpheresFactory.cpp 2010-04-26 13:58:23 +0000
@@ -143,9 +143,6 @@
shared_ptr<Aabb> aabb(new Aabb);
shared_ptr<Sphere> iSphere(new Sphere);
- Quaternionr q;
- q.FromAxisAngle( Vector3r(0,0,1),0);
-
body->isDynamic = true;
physics->velocity = Vector3r(//
@@ -161,7 +158,7 @@
2.0/5.0*physics->mass*r*r,
2.0/5.0*physics->mass*r*r,
2.0/5.0*physics->mass*r*r);
- physics->se3 = Se3r(position,q);
+ physics->se3 = Se3r(position,Quaternionr::Identity());
if (young) physics->young = young;
if (poisson) physics->poisson = poisson;
if (frictionAngle) physics->frictionAngle = frictionAngle;
=== modified file 'pkg/common/Engine/PartialEngine/JumpChangeSe3.cpp'
--- pkg/common/Engine/PartialEngine/JumpChangeSe3.cpp 2010-03-20 12:40:44 +0000
+++ pkg/common/Engine/PartialEngine/JumpChangeSe3.cpp 2010-04-26 13:58:23 +0000
@@ -10,9 +10,9 @@
if(setVelocities){
Real dt=Omega::instance().getTimeStep();
b->state->vel=deltaSe3.position/dt;
- Vector3r axis; Real angle; deltaSe3.orientation.ToAxisAngle(axis,angle); axis.Normalize();
- b->state->angVel=axis*angle/dt;
- LOG_DEBUG("Angular velocity set to "<<axis*angle/dt<<". Axis="<<axis<<", angle="<<angle);
+ AngleAxisr aa(angleAxisFromQuat(deltaSe3.orientation)); aa.axis().normalize();
+ b->state->angVel=aa.axis()*aa.angle()/dt;
+ LOG_DEBUG("Angular velocity set to "<<aa.axis()*aa.angle()/dt<<". Axis="<<aa.axis()<<", angle="<<aa.angle());
}
if(!setVelocities || (setVelocities && !b->isDynamic)){
b->state->pos+=deltaSe3.position;
=== modified file 'pkg/common/Engine/PartialEngine/RotationEngine.cpp'
--- pkg/common/Engine/PartialEngine/RotationEngine.cpp 2010-03-20 12:40:44 +0000
+++ pkg/common/Engine/PartialEngine/RotationEngine.cpp 2010-04-26 13:58:23 +0000
@@ -21,8 +21,7 @@
void SpiralEngine::action(){
Real dt=Omega::instance().getTimeStep();
axis.Normalize();
- Quaternionr q;
- q.FromAxisAngle(axis,angularVelocity*dt);
+ Quaternionr q(AngleAxisr(angularVelocity*dt,axis));
angleTurned+=angularVelocity*dt;
shared_ptr<BodyContainer> bodies = scene->bodies;
FOREACH(body_id_t id,subscribedBodies){
@@ -46,8 +45,7 @@
void RotationEngine::action(){
rotationAxis.Normalize();
- Quaternionr q;
- q.FromAxisAngle(rotationAxis,angularVelocity*scene->dt);
+ Quaternionr q(AngleAxisr(angularVelocity*scene->dt,rotationAxis));
#ifdef YADE_OPENMP
const long size=subscribedBodies.size();
#pragma omp parallel for schedule(static)
=== modified file 'pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp'
--- pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp 2010-04-25 13:18:11 +0000
+++ pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp 2010-04-26 13:58:23 +0000
@@ -81,9 +81,9 @@
if(!b || !b->shape) continue;
glPushMatrix();
const Se3r& se3=b->state->se3;
- Real angle; Vector3r axis; se3.orientation.ToAxisAngle(axis,angle);
+ AngleAxisr aa(angleAxisFromQuat(se3.orientation));
glTranslatef(se3.position[0],se3.position[1],se3.position[2]);
- glRotatef(angle*Mathr::RAD_TO_DEG,axis[0],axis[1],axis[2]);
+ glRotatef(aa.angle()*Mathr::RAD_TO_DEG,aa.axis()[0],aa.axis()[1],aa.axis()[2]);
//if(b->shape->getClassName() != "LineSegment"){ // FIXME: a body needs to say: I am selectable ?!?!
glPushName(b->getId());
shapeDispatcher(b->shape,b->state,wire || b->shape->wire,viewInfo);
@@ -119,9 +119,9 @@
if(!scaleRotations) bodyDisp[id].ori=ori;
else{
Quaternionr relRot=refOri.Conjugate()*ori;
- Vector3r axis; Real angle; relRot.ToAxisAngle(axis,angle);
- angle*=rotScale;
- bodyDisp[id].ori=refOri*Quaternionr(axis,angle);
+ AngleAxisr aa(angleAxisFromQuat(relRot));
+ aa.angle()*=rotScale;
+ bodyDisp[id].ori=refOri*Quaternionr(aa);
}
}
}
=== modified file 'pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp'
--- pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2010-04-25 15:46:26 +0000
+++ pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2010-04-26 13:58:23 +0000
@@ -44,7 +44,7 @@
Quaternionr normal2pt;
Vector3r axis=planeNormal.cross(planePt); axis.Normalize();
Real angle=planePt.norm()/radius;
- normal2pt.FromAxisAngle(axis,angle);
+ normal2pt=Quaternionr(AngleAxisr(angle,axis));
Quaternionr ret; return ret.setFromTwoVectors(Vector3r::UnitX(),normal2pt*planeNormal);
} else {
Quaternionr ret; return ret.setFromTwoVectors(Vector3r::UnitX(),planeNormal);
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp 2010-04-25 15:46:26 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp 2010-04-26 13:58:23 +0000
@@ -56,7 +56,7 @@
//1st imagine we have to rotate shearForce ON the x-y plane as a result of contact point rotation
//axis = phys->prevNormal.Cross(geom->normal); // axis of rotation of contact point
//angle = Mathr::ACos(geom->normal.Dot(phys->prevNormal)); //angle of rotation of contact point
- //q.FromAxisAngle(axis,angle); //quaternion for rotation of contact point
+ //q=Quaternion(AngleAxisr(angleAxis)); //quaternion for rotation of contact point
// shearForce = q*shearForce; //rotate shearForce from previous contact orientation to current contact orientation
//Then imagine we have to find a direction by rotating THROUGH the x-y plane, ABOUT the contact normal (i.e. axis). The direction should be along the resultant angular velocity.
=== modified file 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp'
--- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp 2010-04-25 15:46:26 +0000
+++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp 2010-04-26 13:58:23 +0000
@@ -203,17 +203,15 @@
Vector3r axis = state->angVel;
if (axis!=Vector3r::Zero()) { //If we have an angular velocity, we make a rotation
- Quaternionr q;
Real angle=axis.norm(); axis/=angle;
- q.FromAxisAngle ( axis,angle*dt );
+ Quaternionr q(AngleAxisr(angle*dt,axis));
state->ori = q*state->ori;
}
if(scene->forces.getMoveRotUsed() && scene->forces.getRot(id)!=Vector3r::Zero()) {
- Quaternionr q;
Vector3r r(scene->forces.getRot(id));
Real norm=r.norm(); r/=norm;
- q.FromAxisAngle(r,norm);
+ Quaternionr q(AngleAxisr(norm,r));
state->ori=q*state->ori;
}
state->ori.Normalize();
@@ -237,8 +235,7 @@
if(scene->forces.getMoveRotUsed() && scene->forces.getRot(id)!=Vector3r::Zero()) {
Vector3r r(scene->forces.getRot(id));
Real norm=r.norm(); r/=norm;
- Quaternionr q;
- q.FromAxisAngle(r,norm);
+ Quaternionr q(AngleAxisr(norm,r));
state->ori=q*state->ori;
}
state->ori.Normalize();
@@ -247,9 +244,9 @@
Quaternionr NewtonIntegrator::DotQ(const Vector3r& angVel, const Quaternionr& Q){
// FIXME: uses index access which has different meaning in Wm3 and Eigen
Quaternionr dotQ;
- dotQ[0] = (-Q[1]*angVel[0]-Q[2]*angVel[1]-Q[3]*angVel[2])/2;
- dotQ[1] = ( Q[0]*angVel[0]-Q[3]*angVel[1]+Q[2]*angVel[2])/2;
- dotQ[2] = ( Q[3]*angVel[0]+Q[0]*angVel[1]-Q[1]*angVel[2])/2;
- dotQ[3] = (-Q[2]*angVel[0]+Q[1]*angVel[1]+Q[0]*angVel[2])/2;
+ dotQ.w() = (-Q.x()*angVel[0]-Q.y()*angVel[1]-Q.z()*angVel[2])/2;
+ dotQ.x() = ( Q.w()*angVel[0]-Q.z()*angVel[1]+Q.y()*angVel[2])/2;
+ dotQ.y() = ( Q.z()*angVel[0]+Q.w()*angVel[1]-Q.x()*angVel[2])/2;
+ dotQ.z() = (-Q.y()*angVel[0]+Q.x()*angVel[1]+Q.w()*angVel[2])/2;
return dotQ;
}
=== modified file 'pkg/dem/Engine/GlobalEngine/VTKRecorder.hpp'
--- pkg/dem/Engine/GlobalEngine/VTKRecorder.hpp 2010-03-20 12:40:44 +0000
+++ pkg/dem/Engine/GlobalEngine/VTKRecorder.hpp 2010-04-26 13:58:23 +0000
@@ -6,12 +6,12 @@
public:
enum {REC_SPHERES=0,REC_FACETS,REC_COLORS,REC_CPM,REC_INTR,REC_VELOCITY,REC_IDS,REC_CLUMPIDS,REC_SENTINEL};
virtual void action();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(VTKRecorder,PeriodicEngine,"Engine recording snapshots of simulation into series of *.vtu files, readable by VTK-based postprocessing programs such as Paraview. Both bodies (spheres and facets) and interactions can be recorded, with various vector/scalar quantities that are defined on them.",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(VTKRecorder,PeriodicEngine,"Engine recording snapshots of simulation into series of *.vtu files, readable by VTK-based postprocessing programs such as Paraview. Both bodies (spheres and facets) and interactions can be recorded, with various vector/scalar quantities that are defined on them.\n\n:yref:`PeriodicEngine.initRun` is initialized to ``True`` automatically.",
((bool,compress,false,"Compress output XML files [experimental]."))
((bool,skipFacetIntr,true,"Skip interactions with facets, when saving interactions"))
((bool,skipNondynamic,false,"Skip non-dynamic spheres (but not facets)."))
((string,fileName,"","Base file name; it will be appended with {spheres,intrs,facets}-243100.vtu depending on active recorders and step number (243100 in this case). It can contain slashes, but the directory must exist already."))
- ((vector<string>,recorders,,"List of active recorders (as strings). Acceptable recorders are spheres, velocity, facets, colors, cpm, intr, ids, clumpids.")),
+ ((vector<string>,recorders,,"List of active recorders (as strings). Accepted recorders are:\n\n``spheres``\n\tSaves positions and radii (``radii``) of :yref:`spherical<Sphere>` particles.\n``ids``\n\tSaves id's (field ``IDS``) of spheres; active only if ``spheres`` is active.\n``clumpids``\n\tSaves id's of clumps to which each sphere belongs (field ``clumpIDS``); active only if ``spheres`` is active.\n``colors``\n\tSaves colors of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``colors`` spheres and ``Colors`` for facets); only active if ``spheres`` or ``facets`` is activated.\n``velocity``\n\tSaves linear and angular velocities of spherical particles (fields ``velocity`` and ``angVel`` respectively``); only effective with ``spheres``.\n``facets``\n\tSave :yref:`facets<Facet>` positions (vertices).\n``cpm``\n\tSaves data pertaining to the :yref:`concrete model<Law2_Dem3DofGeom_CpmPhys_Cpm>`: ``cpmDamage`` (normalized residual strength averaged on particle), ``cpmSigma`` (stress on particle, normal components), ``cpmTau`` (shear components of stress on particle), ``cpmSigmaM`` (mean stress around particle); ``intr`` is activated automatically by ``cpm``.\n``intr``\n\tWhen ``cpm`` is used, it saves magnitude of normal (``forceN``) and shear (``absForceT``) forces.\n\n\tWithout ``cpm``, saves [TODO]\n\n")),
/*ctor*/ initRun=true;
);
DECLARE_LOGGER;
=== modified file 'pkg/dem/Engine/PartialEngine/Disp2DPropLoadEngine.cpp'
--- pkg/dem/Engine/PartialEngine/Disp2DPropLoadEngine.cpp 2010-04-25 15:46:26 +0000
+++ pkg/dem/Engine/PartialEngine/Disp2DPropLoadEngine.cpp 2010-04-26 13:58:23 +0000
@@ -124,8 +124,7 @@
dalpha = Mathr::ATan( (A - Mathr::Tan(alpha))/(1.0 + A * Mathr::Tan(alpha)));
}
- Quaternionr qcorr;
- qcorr.FromAxisAngle(Vector3r(0,0,1),dalpha);
+ Quaternionr qcorr(AngleAxisr(dalpha,Vector3r::UnitZ()));
if(LOG)
cout << "Quaternion associe a la rotation incrementale : " << qcorr.w() << " " << qcorr.x() << " " << qcorr.y() << " " << qcorr.z() << endl;
=== modified file 'pkg/dem/Engine/PartialEngine/KinemCNDEngine.cpp'
--- pkg/dem/Engine/PartialEngine/KinemCNDEngine.cpp 2010-04-03 16:40:33 +0000
+++ pkg/dem/Engine/PartialEngine/KinemCNDEngine.cpp 2010-04-26 13:58:23 +0000
@@ -81,8 +81,8 @@
dalpha = Mathr::ATan( (A - Mathr::Tan(alpha))/(1.0 + A * Mathr::Tan(alpha)));
}
- Quaternionr qcorr;
- qcorr.FromAxisAngle(Vector3r(0,0,1),dalpha);
+ Quaternionr qcorr(AngleAxisr(dalpha,Vector3r::UnitZ()));
+
// On applique la rotation en changeant l'orientation des plaques, leurs vang et en affectant donc alpha
leftbox->state->ori = qcorr*leftbox->state->ori;
=== modified file 'pkg/dem/Engine/PartialEngine/KinemCNLEngine.cpp'
--- pkg/dem/Engine/PartialEngine/KinemCNLEngine.cpp 2010-04-03 16:40:33 +0000
+++ pkg/dem/Engine/PartialEngine/KinemCNLEngine.cpp 2010-04-26 13:58:23 +0000
@@ -107,9 +107,8 @@
Real A = (Ysup_mod - Ylat_mod) * 2.0*Mathr::Tan(alpha) / (2.0*(Ysup - Ylat) + dx*Mathr::Tan(alpha) );
dalpha = Mathr::ATan( (A - Mathr::Tan(alpha))/(1.0 + A * Mathr::Tan(alpha)));
}
-
- Quaternionr qcorr;
- qcorr.FromAxisAngle(Vector3r(0,0,1),dalpha);
+
+ Quaternionr qcorr(AngleAxisr(dalpha,Vector3r::UnitZ()));
// On applique la rotation en changeant l'orientation des plaques, leurs vang et en affectant donc alpha
leftbox->state->ori = qcorr*leftbox->state->ori;
=== modified file 'pkg/dem/Engine/PartialEngine/KinemCNSEngine.cpp'
--- pkg/dem/Engine/PartialEngine/KinemCNSEngine.cpp 2010-04-03 16:40:33 +0000
+++ pkg/dem/Engine/PartialEngine/KinemCNSEngine.cpp 2010-04-26 13:58:23 +0000
@@ -92,8 +92,7 @@
dalpha = Mathr::ATan( (A - Mathr::Tan(alpha))/(1.0 + A * Mathr::Tan(alpha)));
}
- Quaternionr qcorr;
- qcorr.FromAxisAngle(Vector3r(0,0,1),dalpha);
+ Quaternionr qcorr(AngleAxisr(dalpha,Vector3r::UnitZ()));
if(LOG)
cout << "Quaternion associe a la rotation incrementale : " << qcorr.W() << " " << qcorr.X() << " " << qcorr.Y() << " " << qcorr.Z() << endl;
// On applique la rotation en changeant l'orientation des plaques, leurs vang et en affectant donc alpha
=== modified file 'pkg/dem/Engine/PartialEngine/KinemCTDEngine.cpp'
--- pkg/dem/Engine/PartialEngine/KinemCTDEngine.cpp 2010-04-25 15:46:26 +0000
+++ pkg/dem/Engine/PartialEngine/KinemCTDEngine.cpp 2010-04-26 13:58:23 +0000
@@ -119,8 +119,7 @@
dalpha = Mathr::ATan( (A - Mathr::Tan(alpha))/(1.0 + A * Mathr::Tan(alpha)));
}
- Quaternionr qcorr;
- qcorr.FromAxisAngle(Vector3r(0,0,1),dalpha);
+ Quaternionr qcorr(AngleAxisr(dalpha,Vector3r::UnitZ()));
// On applique la rotation en changeant l'orientation des plaques et leurs vang
leftbox->state->ori = qcorr*leftbox->state->ori;
=== modified file 'pkg/dem/PreProcessor/CohesiveTriaxialTest.cpp'
--- pkg/dem/PreProcessor/CohesiveTriaxialTest.cpp 2010-04-25 15:46:26 +0000
+++ pkg/dem/PreProcessor/CohesiveTriaxialTest.cpp 2010-04-26 13:58:23 +0000
@@ -285,7 +285,6 @@
Quaternionr q(Mathr::SymmetricRandom(),Mathr::SymmetricRandom(),Mathr::SymmetricRandom(),Mathr::SymmetricRandom());
q.Normalize();
-// q.FromAxisAngle( Vector3r(0,0,1),0);
body->isDynamic = dynamic;
@@ -328,9 +327,6 @@
shared_ptr<Box> iBox(new Box);
- Quaternionr q;
- q.FromAxisAngle( Vector3r(0,0,1),0);
-
body->isDynamic = false;
body->state->angVel = Vector3r(0,0,0);
@@ -344,7 +340,7 @@
);
// physics->mass = 0;
// physics->inertia = Vector3r(0,0,0);
- body->state->se3 = Se3r(position,q);
+ body->state->se3 = Se3r(position,Quaternionr::Identity());
physics->young = boxYoungModulus;
physics->poisson = boxKsDivKn;
@@ -451,9 +447,6 @@
void CohesiveTriaxialTest::positionRootBody(shared_ptr<Scene>& rootBody)
{
- Quaternionr q;
- q.FromAxisAngle( Vector3r(0,0,1),0);
-
shared_ptr<Aabb> aabb(new Aabb);
aabb->diffuseColor = Vector3r(0,0,1);
}
=== modified file 'pkg/dem/PreProcessor/DirectShearCis.cpp'
--- pkg/dem/PreProcessor/DirectShearCis.cpp 2010-04-25 13:18:11 +0000
+++ pkg/dem/PreProcessor/DirectShearCis.cpp 2010-04-26 13:58:23 +0000
@@ -175,16 +175,13 @@
shared_ptr<SphereModel> gSphere(new SphereModel);
shared_ptr<Sphere> iSphere(new Sphere);
- Quaternionr q;
- q.FromAxisAngle( Vector3r(0,0,1),0);
-
body->isDynamic = true;
physics->angularVelocity = Vector3r(0,0,0);
physics->velocity = Vector3r(0,0,0);
physics->mass = 4.0/3.0*Mathr::PI*radius*radius*radius*density;
physics->inertia = Vector3r(2.0/5.0*physics->mass*radius*radius,2.0/5.0*physics->mass*radius*radius,2.0/5.0*physics->mass*radius*radius);
- physics->se3 = Se3r(position,q);
+ physics->se3 = Se3r(position,Quaternionr::Identity());
physics->young = sphereYoungModulus;
physics->poisson = spherePoissonRatio;
physics->frictionAngle = sphereFrictionDeg * Mathr::PI/180.0;
@@ -217,9 +214,6 @@
shared_ptr<Box> iBox(new Box);
- Quaternionr q;
- q.FromAxisAngle( Vector3r(0,0,1),0);
-
body->isDynamic = false;
physics->angularVelocity = Vector3r(0,0,0);
@@ -232,7 +226,7 @@
);
//physics->mass = 0;
//physics->inertia = Vector3r(0,0,0);
- physics->se3 = Se3r(position,q);
+ physics->se3 = Se3r(position,Quaternion::Identity());
physics->young = boxYoungModulus;
physics->poisson = boxPoissonRatio;
physics->frictionAngle = 0.0; //default value, modified after for w2 and w4 to have good values of phi(sphere-walls)
@@ -338,11 +332,8 @@
{
rootBody->isDynamic = false;
- Quaternionr q;
- q.FromAxisAngle( Vector3r(0,0,1),0);
-
shared_ptr<ParticleParameters> physics(new ParticleParameters); // FIXME : fix indexable class PhysicalParameters
- physics->se3 = Se3r(Vector3r(0,0,0),q);
+ physics->se3 = Se3r(Vector3r(0,0,0),Quaternionr::Identity());
physics->mass = 0;
physics->velocity = Vector3r(0,0,0);
physics->acceleration = Vector3r::Zero();
=== modified file 'pkg/dem/PreProcessor/STLImporterTest.cpp'
--- pkg/dem/PreProcessor/STLImporterTest.cpp 2010-04-25 13:18:11 +0000
+++ pkg/dem/PreProcessor/STLImporterTest.cpp 2010-04-26 13:58:23 +0000
@@ -150,9 +150,6 @@
shared_ptr<Aabb> aabb(new Aabb);
shared_ptr<Sphere> iSphere(new Sphere);
- Quaternionr q;
- q.FromAxisAngle( Vector3r(0,0,1),0);
-
Vector3r position = Vector3r(i,j+spheresHeight,k)*(2*maxRadius*1.1)
- Vector3r( nbSpheres[0]/2*(2*maxRadius*1.1) , 0 , nbSpheres[2]/2*(2*maxRadius*1.1) )
+ Vector3r( Mathr::SymmetricRandom()*disorder[0]
@@ -167,7 +164,7 @@
physics->velocity = Vector3r(0,0,0);
physics->mass = 4.0/3.0*Mathr::PI*radius*radius*radius*density;
physics->inertia = Vector3r(2.0/5.0*physics->mass*radius*radius,2.0/5.0*physics->mass*radius*radius,2.0/5.0*physics->mass*radius*radius); //
- physics->se3 = Se3r(position,q);
+ physics->se3 = Se3r(position,Quaternionr::Identity());
physics->young = sphereYoungModulus;
physics->poisson = spherePoissonRatio;
physics->frictionAngle = sphereFrictionDeg * Mathr::PI/180.0;
@@ -256,11 +253,8 @@
{
rootBody->isDynamic = false;
- Quaternionr q;
- q.FromAxisAngle( Vector3r(0,0,1),0);
-
shared_ptr<ParticleParameters> physics(new ParticleParameters); // FIXME : fix indexable class PhysicalParameters
- physics->se3 = Se3r(Vector3r(0,0,0),q);
+ physics->se3 = Se3r(Vector3r(0,0,0),Quaternionr::Identity());
physics->mass = 0;
physics->velocity = Vector3r(0,0,0);
physics->acceleration = Vector3r::Zero();
=== modified file 'pkg/dem/PreProcessor/SimpleShear.cpp'
--- pkg/dem/PreProcessor/SimpleShear.cpp 2010-04-25 15:46:26 +0000
+++ pkg/dem/PreProcessor/SimpleShear.cpp 2010-04-26 13:58:23 +0000
@@ -164,12 +164,9 @@
// shared_ptr<SphereModel> gSphere(new SphereModel);
shared_ptr<Sphere> iSphere(new Sphere);
- Quaternionr q; // to define the initial orientation of the sphere
- q.FromAxisAngle( Vector3r(0,0,1),0);
-
body->isDynamic = true;
body->state->pos =position;
- body->state->ori =q;
+ body->state->ori =Quaternionr::Identity();
body->state->vel =Vector3r(0,0,0);
body->state->angVel =Vector3r(0,0,0);
@@ -209,18 +206,13 @@
// shared_ptr<BoxModel> gBox(new BoxModel);
shared_ptr<Box> iBox(new Box);
-
-
- Quaternionr q;
- q.FromAxisAngle( Vector3r(0,0,1),0);
-
body->isDynamic = false;
body->state->angVel = Vector3r(0,0,0);
body->state->vel = Vector3r(0,0,0);
// NB : mass and inertia not defined because not used, since Box are not dynamics
body->state->pos = position;
- body->state->ori = q;
+ body->state->ori = Quaternionr::Identity();
mat->young = boxYoungModulus;
mat->poisson = boxPoissonRatio;
=== modified file 'pkg/dem/PreProcessor/TriaxialTest.cpp'
--- pkg/dem/PreProcessor/TriaxialTest.cpp 2010-04-25 15:46:26 +0000
+++ pkg/dem/PreProcessor/TriaxialTest.cpp 2010-04-26 13:58:23 +0000
@@ -220,8 +220,6 @@
shared_ptr<Aabb> aabb(new Aabb);
shared_ptr<Sphere> iSphere(new Sphere);
- Quaternionr q;
- q.FromAxisAngle( Vector3r(0,0,1),0);
body->isDynamic = dynamic;
body->state->mass = 4.0/3.0*Mathr::PI*radius*radius*radius*density;
body->state->inertia = Vector3r( 2.0/5.0*body->state->mass*radius*radius,
=== modified file 'pkg/dem/PreProcessor/TriaxialTestWater.cpp'
--- pkg/dem/PreProcessor/TriaxialTestWater.cpp 2010-04-25 15:46:26 +0000
+++ pkg/dem/PreProcessor/TriaxialTestWater.cpp 2010-04-26 13:58:23 +0000
@@ -346,9 +346,6 @@
shared_ptr<Aabb> aabb(new Aabb);
shared_ptr<Box> iBox(new Box);
- Quaternionr q;
- q.FromAxisAngle( Vector3r(0,0,1),0);
-
body->isDynamic = false;
body->state->angVel = Vector3r(0,0,0);
@@ -362,7 +359,7 @@
);
// physics->mass = 0;
// physics->inertia = Vector3r(0,0,0);
- body->state->se3 = Se3r(position,q);
+ body->state->se3 = Se3r(position,Quaternionr::Identity());
physics->young = boxYoungModulus;
physics->poisson = boxKsDivKn;
=== modified file 'py/_eudoxos.cpp'
--- py/_eudoxos.cpp 2010-04-25 13:18:11 +0000
+++ py/_eudoxos.cpp 2010-04-26 13:58:23 +0000
@@ -187,7 +187,7 @@
CpmPhys* phys=YADE_CAST<CpmPhys*>(fi.i->interactionPhysics.get());
// transformation matrix, to rotate to the plane
Vector3r ax(Vector3r::Zero()); ax[axis]=1.;
- Quaternionr q; q.FromAxisAngle(ax,fi.theta); q=q.Conjugate();
+ Quaternionr q(AngleAxisr(fi.theta,ax)); q=q.Conjugate();
Matrix3r TT; q.ToRotationMatrix(TT);
//
Real d=(geom->se31.position-geom->se32.position).Length(); // current contact length
=== modified file 'py/yadeWrapper/customConverters.cpp'
--- py/yadeWrapper/customConverters.cpp 2010-03-16 15:34:25 +0000
+++ py/yadeWrapper/customConverters.cpp 2010-04-26 13:58:23 +0000
@@ -88,7 +88,7 @@
se3->position=Vector3r(extract<Real>(PySequence_GetItem(obj_ptr,0)),extract<Real>(PySequence_GetItem(obj_ptr,1)),extract<Real>(PySequence_GetItem(obj_ptr,2)));
Vector3r axis=Vector3r(extract<Real>(PySequence_GetItem(obj_ptr,3)),extract<Real>(PySequence_GetItem(obj_ptr,4)),extract<Real>(PySequence_GetItem(obj_ptr,5)));
Real angle=extract<Real>(PySequence_GetItem(obj_ptr,6));
- se3->orientation.FromAxisAngle(axis,angle);
+ se3->orientation=Quaternionr(AngleAxisr(angle,axis));
} else throw std::logic_error(__FILE__ ": First, the sequence size for Se3r object was 2 or 7, but now is not? (programming error, please report!");
data->convertible=storage;
}
=== modified file 'py/yadeWrapper/yadeWrapper.cpp'
--- py/yadeWrapper/yadeWrapper.cpp 2010-04-25 11:58:30 +0000
+++ py/yadeWrapper/yadeWrapper.cpp 2010-04-26 13:58:23 +0000
@@ -295,7 +295,7 @@
double simulationTime(){return OMEGA.getSimulationTime();}
double realTime(){ return OMEGA.getComputationTime(); }
double dt_get(){return OMEGA.getTimeStep();}
- void dt_set(double dt){OMEGA.skipTimeStepper(true); OMEGA.setTimeStep(dt); LOG_WARN("Setting timestep "<<dt<<", "<<dt_get()); }
+ void dt_set(double dt){OMEGA.skipTimeStepper(true); OMEGA.setTimeStep(dt); }
long stopAtIter_get(){return OMEGA.getScene()->stopAtIteration; }
void stopAtIter_set(long s){OMEGA.getScene()->stopAtIteration=s; }
=== modified file 'scripts/test-spiral.py'
--- scripts/test-spiral.py 2010-03-22 17:39:33 +0000
+++ scripts/test-spiral.py 2010-04-26 13:58:23 +0000
@@ -1,9 +1,9 @@
# script for testing InterpolatingSpiralEngine: sphere going in a sphere-like motion around bar
O.bodies.append([utils.box([0,0,0],[.005,.005,1],dynamic=False),utils.sphere([0,.1,-1],.04,dynamic=False)])
O.engines=[
- InterpolatingSpiralEngine(subscribedBodies=[1],times=[10,20,30,40,50,60,70,80,90,100],angularVelocities=[1,2,3,4,5,3,1,-1,-3,0],axis=[0,0,1],axisPt=[0,0,0],period=-1,slope=.003,label='spiral'),
+ InterpolatingSpiralEngine(subscribedBodies=[1],times=[10,20,30,40,50,60,70,80,90,100],angularVelocities=[1,2,3,4,5,3,1,-1,-3,0],axis=[0,0,1],axisPt=[0,0,0],wrap=True,slope=.003,label='spiral'),
]
-O.dt=5e-4
+O.dt=4e-6
O.saveTmp('initial')
from yade import qt
qt.Controller()
Follow ups
