yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #04824
[Branch ~yade-dev/yade/trunk] Rev 2286: 1. Replace wm3-compatibility functions (component{Min, Max}Vector, angleAxisFromQuat, diag{Mult, Di...
------------------------------------------------------------
revno: 2286
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Wed 2010-06-09 00:25:00 +0200
message:
1. Replace wm3-compatibility functions (component{Min,Max}Vector, angleAxisFromQuat, diag{Mult,Div}) with their eigen equivalents.
modified:
core/Scene.cpp
core/State.hpp
gui/qt3/GLViewer.cpp
lib/base/Math.hpp
lib/serialization/KnownFundamentalsHandler.tpp
pkg/common/Engine/Functor/Bo1_Facet_Aabb.cpp
pkg/common/Engine/PartialEngine/StepDisplacer.cpp
pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp
pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp
pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp
pkg/dem/DataClass/SpherePack.hpp
pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp
pkg/dem/Engine/GlobalEngine/FacetTopologyAnalyzer.cpp
pkg/dem/Engine/GlobalEngine/MicroMacroAnalyser.cpp
pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp
pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp
pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp
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/meta/ConcretePM.cpp
pkg/dem/meta/Shop.cpp
py/_eudoxos.cpp
py/_utils.cpp
py/mathWrap/miniEigen.cpp
py/pack/_packObb.cpp
py/pack/_packPredicates.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/Scene.cpp'
--- core/Scene.cpp 2010-05-13 20:19:38 +0000
+++ core/Scene.cpp 2010-06-08 22:25:00 +0000
@@ -159,8 +159,8 @@
if(!isinf(b->bound->min[i])) mn[i]=min(mn[i],b->bound->min[i]);
}
} else {
- mx=componentMaxVector(mx,b->state->pos);
- mn=componentMinVector(mn,b->state->pos);
+ mx=mx.cwise().max(b->state->pos);
+ mn=mn.cwise().min(b->state->pos);
}
}
bound->min=mn; bound->max=mx;
=== modified file 'core/State.hpp'
--- core/State.hpp 2010-06-03 17:49:44 +0000
+++ core/State.hpp 2010-06-08 22:25:00 +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; AngleAxisr aa(angleAxisFromQuat(relRot)); return aa.axis()*aa.angle(); }
+ Vector3r rot() const { Quaternionr relRot=refOri.conjugate()*ori; AngleAxisr aa(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-06-08 21:22:07 +0000
+++ gui/qt3/GLViewer.cpp 2010-06-08 22:25:00 +0000
@@ -402,8 +402,8 @@
Real inf=std::numeric_limits<Real>::infinity();
min=Vector3r(inf,inf,inf); max=Vector3r(-inf,-inf,-inf);
FOREACH(const shared_ptr<Body>& b, *rb->bodies){
- max=componentMaxVector(max,b->state->pos);
- min=componentMinVector(min,b->state->pos);
+ max=max.cwise().max(b->state->pos);
+ min=min.cwise().min(b->state->pos);
}
} else {LOG_DEBUG("Using rootBody's Aabb");}
} else {
@@ -606,7 +606,7 @@
if(!renderer->clipPlaneActive[planeId] && planeId!=manipulatedClipPlane) continue;
glPushMatrix();
const Se3r& se3=renderer->clipPlaneSe3[planeId];
- AngleAxisr aa(angleAxisFromQuat(se3.orientation));
+ AngleAxisr aa(se3.orientation);
glTranslatef(se3.position[0],se3.position[1],se3.position[2]);
glRotated(aa.angle()*Mathr::RAD_TO_DEG,aa.axis()[0],aa.axis()[1],aa.axis()[2]);
Real cff=1;
=== modified file 'lib/base/Math.hpp'
--- lib/base/Math.hpp 2010-06-07 17:09:08 +0000
+++ lib/base/Math.hpp 2010-06-08 22:25:00 +0000
@@ -70,8 +70,6 @@
//template<class Scalar> MATRIX3_TEMPLATE(Scalar) operator*(Scalar s, const MATRIX3_TEMPLATE(Scalar)& m) { return m*s; }
//template<class Scalar> Quaternion<Scalar> operator*(Scalar s, const Quaternion<Scalar>& q) { return q*s; }
template<typename Scalar> void matrixEigenDecomposition(const MATRIX3_TEMPLATE(Scalar) m, MATRIX3_TEMPLATE(Scalar)& mRot, MATRIX3_TEMPLATE(Scalar)& mDiag){ Eigen::SelfAdjointEigenSolver<MATRIX3_TEMPLATE(Scalar)> a(m); mRot=a.eigenvectors(); mDiag=a.eigenvalues().asDiagonal(); }
-//__attribute__((warning("Replace this function with direct AngleAxis constructor from Quaternion")))
-template<typename Scalar> AngleAxis<Scalar> angleAxisFromQuat(const Quaternion<Scalar>& q){ return AngleAxis<Scalar>(q); }
// http://eigen.tuxfamily.org/dox/TutorialGeometry.html
template<typename Scalar> MATRIX3_TEMPLATE(Scalar) matrixFromEulerAnglesXYZ(Scalar x, Scalar y, Scalar z){ MATRIX3_TEMPLATE(Scalar) m; m=AngleAxis<Scalar>(x,VECTOR3_TEMPLATE(Scalar)::UnitX())*AngleAxis<Scalar>(y,VECTOR3_TEMPLATE(Scalar)::UnitY())*AngleAxis<Scalar>(z,VECTOR3_TEMPLATE(Scalar)::UnitZ()); return m;}
template<typename Scalar> bool operator==(const Quaternion<Scalar>& u, const Quaternion<Scalar>& v){ return u.x()==v.x() && u.y()==v.y() && u.z()==v.z() && u.w()==v.w(); }
@@ -118,32 +116,6 @@
* TODO
*/
-// replace by outer product of 2 vectors: v1*v2.transpose();
-template<typename Scalar>
-MATRIX3_TEMPLATE(Scalar) makeTensorProduct (const VECTOR3_TEMPLATE(Scalar)& rkU, const VECTOR3_TEMPLATE(Scalar)& rkV)
-{
- MATRIX3_TEMPLATE(Scalar) ret;
- ret(0,0) = rkU[0]*rkV[0];
- ret(0,1) = rkU[0]*rkV[1];
- ret(0,2) = rkU[0]*rkV[2];
- ret(1,0) = rkU[1]*rkV[0];
- ret(1,1) = rkU[1]*rkV[1];
- ret(1,2) = rkU[1]*rkV[2];
- ret(2,0) = rkU[2]*rkV[0];
- ret(2,1) = rkU[2]*rkV[1];
- ret(2,2) = rkU[2]*rkV[2];
- return ret;
-}
-// eigen2: a.cwise().min(b) / max(b)
-// eigen3: a.array().min(b) / max(b)
-template<typename Scalar> VECTOR2_TEMPLATE(Scalar) componentMaxVector(const VECTOR2_TEMPLATE(Scalar)& a, const VECTOR2_TEMPLATE(Scalar)& b){ return VECTOR2_TEMPLATE(Scalar)(std::max(a.x(),b.x()),std::max(a.y(),b.y()));}
-template<typename Scalar> VECTOR2_TEMPLATE(Scalar) componentMinVector(const VECTOR2_TEMPLATE(Scalar)& a, const VECTOR2_TEMPLATE(Scalar)& b){ return VECTOR2_TEMPLATE(Scalar)(std::min(a.x(),b.x()),std::min(a.y(),b.y()));}
-template<typename Scalar> VECTOR3_TEMPLATE(Scalar) componentMaxVector(const VECTOR3_TEMPLATE(Scalar)& a, const VECTOR3_TEMPLATE(Scalar)& b){ return VECTOR3_TEMPLATE(Scalar)(std::max(a.x(),b.x()),std::max(a.y(),b.y()),std::max(a.z(),b.z()));}
-template<typename Scalar> VECTOR3_TEMPLATE(Scalar) componentMinVector(const VECTOR3_TEMPLATE(Scalar)& a, const VECTOR3_TEMPLATE(Scalar)& b){ return VECTOR3_TEMPLATE(Scalar)(std::min(a.x(),b.x()),std::min(a.y(),b.y()),std::min(a.z(),b.z()));}
-// eigen2: v1.cwise()*v2;
-// eigen3: v1.array()*v2.array()
-template<typename Scalar> VECTOR3_TEMPLATE(Scalar) diagMult(const VECTOR3_TEMPLATE(Scalar)& a, const VECTOR3_TEMPLATE(Scalar)& b){return VECTOR3_TEMPLATE(Scalar)(a.x()*b.x(),a.y()*b.y(),a.z()*b.z());}
-template<typename Scalar> VECTOR3_TEMPLATE(Scalar) diagDiv(const VECTOR3_TEMPLATE(Scalar)& a, const VECTOR3_TEMPLATE(Scalar)& b){return VECTOR3_TEMPLATE(Scalar)(a.x()/b.x(),a.y()/b.y(),a.z()/b.z());}
// eigen: m << m00,m01,m02,m10,m11,m12,m20,m21,m22;
template<typename Scalar> MATRIX3_TEMPLATE(Scalar) matrixFromElements(Scalar m00, Scalar m01, Scalar m02, Scalar m10, Scalar m11, Scalar m12, Scalar m20, Scalar m21, Scalar m22){ MATRIX3_TEMPLATE(Scalar) m; m(0,0)=m00; m(0,1)=m01; m(0,2)=m02; m(1,0)=m10; m(1,1)=m11; m(1,2)=m12; m(2,0)=m20; m(2,1)=m21; m(2,2)=m22; return m; }
=== modified file 'lib/serialization/KnownFundamentalsHandler.tpp'
--- lib/serialization/KnownFundamentalsHandler.tpp 2010-05-04 13:56:05 +0000
+++ lib/serialization/KnownFundamentalsHandler.tpp 2010-06-08 22:25:00 +0000
@@ -354,7 +354,7 @@
string * tmpStr = any_cast<string*>(a);
Quaternion<RealType> * tmp = any_cast<Quaternion<RealType>*>(ac.getAddress());
- AngleAxis<RealType> aa(angleAxisFromQuat(*tmp));
+ AngleAxis<RealType> aa(*tmp);
aa.axis().normalize();
*tmpStr = IOFormatManager::getCustomFundamentalOpeningBracket() +
@@ -457,7 +457,7 @@
VECTOR3_TEMPLATE(RealType) position;
- AngleAxis<RealType> aa(angleAxisFromQuat(tmp->orientation));
+ AngleAxis<RealType> aa(tmp->orientation);
aa.axis().normalize();
position = tmp->position;
=== modified file 'pkg/common/Engine/Functor/Bo1_Facet_Aabb.cpp'
--- pkg/common/Engine/Functor/Bo1_Facet_Aabb.cpp 2010-05-04 13:56:05 +0000
+++ pkg/common/Engine/Functor/Bo1_Facet_Aabb.cpp 2010-06-08 22:25:00 +0000
@@ -25,16 +25,16 @@
for (int i=1;i<3;++i)
{
Vector3r v = O + facetAxisT * vertices[i];
- aabb->min = componentMinVector( aabb->min, v);
- aabb->max = componentMaxVector( aabb->max, v);
+ aabb->min = aabb->min.cwise().min(v);
+ aabb->max = aabb->max.cwise().max(v);
}
} else {
Real inf=std::numeric_limits<Real>::infinity();
aabb->min=Vector3r(inf,inf,inf); aabb->max=Vector3r(-inf,-inf,-inf);
for(int i=0; i<3; i++){
Vector3r v=scene->cell->unshearPt(O+facetAxisT*vertices[i]);
- aabb->min=componentMinVector(aabb->min,v);
- aabb->max=componentMaxVector(aabb->max,v);
+ aabb->min=aabb->min.cwise().min(v);
+ aabb->max=aabb->min.cwise().max(v);
}
}
}
=== modified file 'pkg/common/Engine/PartialEngine/StepDisplacer.cpp'
--- pkg/common/Engine/PartialEngine/StepDisplacer.cpp 2010-06-08 21:22:07 +0000
+++ pkg/common/Engine/PartialEngine/StepDisplacer.cpp 2010-06-08 22:25:00 +0000
@@ -11,7 +11,7 @@
if(setVelocities){
const Real& dt=scene->dt;
b->state->vel=deltaSe3.position/dt;
- AngleAxisr aa(angleAxisFromQuat(deltaSe3.orientation)); aa.axis().normalize();
+ AngleAxisr aa(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());
}
=== modified file 'pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp'
--- pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp 2010-05-04 13:56:05 +0000
+++ pkg/common/RenderingEngine/OpenGLRenderingEngine.cpp 2010-06-08 22:25:00 +0000
@@ -81,7 +81,7 @@
if(!b || !b->shape) continue;
glPushMatrix();
const Se3r& se3=b->state->se3;
- AngleAxisr aa(angleAxisFromQuat(se3.orientation));
+ AngleAxisr aa(se3.orientation);
glTranslatef(se3.position[0],se3.position[1],se3.position[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 ?!?!
@@ -115,11 +115,11 @@
if(!(scaleDisplacements||scaleRotations||scene->isPeriodic)){ bodyDisp[id].pos=pos; bodyDisp[id].ori=ori; continue; }
// apply scaling
bodyDisp[id].pos=cellPos; // point of reference (inside the cell for periodic)
- if(scaleDisplacements) bodyDisp[id].pos+=diagMult(dispScale,Vector3r(pos-refPos)); // add scaled translation to the point of reference
+ if(scaleDisplacements) bodyDisp[id].pos+=dispScale.cwise()*Vector3r(pos-refPos); // add scaled translation to the point of reference
if(!scaleRotations) bodyDisp[id].ori=ori;
else{
Quaternionr relRot=refOri.conjugate()*ori;
- AngleAxisr aa(angleAxisFromQuat(relRot));
+ AngleAxisr aa(relRot);
aa.angle()*=rotScale;
bodyDisp[id].ori=refOri*Quaternionr(aa);
}
@@ -132,7 +132,7 @@
glColor3v(Vector3r(1,1,0));
glPushMatrix();
Vector3r size=scene->cell->getSize();
- if(dispScale!=Vector3r::Ones()) size+=diagMult(dispScale,Vector3r(size-scene->cell->refSize));
+ if(dispScale!=Vector3r::Ones()) size+=dispScale.cwise()*Vector3r(size-scene->cell->refSize);
glTranslatev(scene->cell->shearPt(.5*size)); // shear center (moves when sheared)
glMultMatrixd(scene->cell->getGlShearTrsfMatrix());
glScalev(size);
@@ -327,7 +327,7 @@
Quaternionr ori=bodyDisp[b->getId()].ori;
if(!b->shape || !((b->getGroupMask()&mask) || b->getGroupMask()==0)) continue;
glPushMatrix();
- AngleAxisr aa(angleAxisFromQuat(ori));
+ AngleAxisr aa(ori);
glTranslatef(pos[0],pos[1],pos[2]);
glRotatef(aa.angle()*Mathr::RAD_TO_DEG,aa.axis()[0],aa.axis()[1],aa.axis()[2]);
if(current_selection==b->getId() || b->shape->highlight){
=== modified file 'pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp'
--- pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp 2010-05-28 09:35:32 +0000
+++ pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_FacetSphere.cpp 2010-06-08 22:25:00 +0000
@@ -130,7 +130,7 @@
#ifdef FACET_TOPO
if(noVertexContact && facet->edgeAdjIds[edgeMax]!=Body::ID_NONE){
// find angle between our normal and the facet's normal (still local coords)
- Quaternionr q; q.Align(facet->nf,normal); AngleAxisr aa(angleAxisFromQuat(q));
+ Quaternionr q; q.Align(facet->nf,normal); AngleAxisr aa(q);
assert(aa.angle()>=0 && aa.angle()<=Mathr::PI);
if(edgeNormals[edgeMax].Dot(aa.axis())<0) aa.angle()*=-1.;
bool negFace=normal.Dot(facet->nf)<0; // contact in on the negative facet's face
=== modified file 'pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp'
--- pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2010-05-28 09:35:32 +0000
+++ pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2010-06-08 22:25:00 +0000
@@ -24,7 +24,7 @@
*/
Vector3r Dem3DofGeom_SphereSphere::unrollSpherePtToPlane(const Quaternionr& fromXtoPtOri, const Real& radius, const Vector3r& planeNormal){
Quaternionr normal2pt; normal2pt.setFromTwoVectors(planeNormal,fromXtoPtOri*Vector3r::UnitX());
- AngleAxisr aa(angleAxisFromQuat(normal2pt));
+ AngleAxisr aa(normal2pt);
return (aa.angle()*radius) /* length */ *(aa.axis().cross(planeNormal)) /* direction: both are unit vectors */;
}
@@ -114,7 +114,7 @@
// FIXME: this is not correct, as it assumes normal will not change (?)
Quaternionr relOri12=ori1.conjugate()*ori2;
Quaternionr oriDiff=initRelOri12.conjugate()*relOri12;
- AngleAxisr aa(angleAxisFromQuat(oriDiff));
+ AngleAxisr aa(oriDiff);
if(aa.angle()>Mathr::PI)aa.angle()-=Mathr::TWO_PI;
// cerr<<axis<<";"<<angle<<";"<<ori1<<";"<<ori2<<";"<<oriDiff<<endl;
return aa.angle()*aa.axis();
=== modified file 'pkg/dem/DataClass/SpherePack.hpp'
--- pkg/dem/DataClass/SpherePack.hpp 2010-05-04 13:56:05 +0000
+++ pkg/dem/DataClass/SpherePack.hpp 2010-06-08 22:25:00 +0000
@@ -62,7 +62,7 @@
python::tuple aabb_py() const { Vector3r mn,mx; aabb(mn,mx); return python::make_tuple(mn,mx); }
void aabb(Vector3r& mn, Vector3r& mx) const {
Real inf=std::numeric_limits<Real>::infinity(); mn=Vector3r(inf,inf,inf); mx=Vector3r(-inf,-inf,-inf);
- FOREACH(const Sph& s, pack){ Vector3r r(s.r,s.r,s.r); mn=componentMinVector(mn,Vector3r(s.c-r)); mx=componentMaxVector(mx,Vector3r(s.c+r));}
+ FOREACH(const Sph& s, pack){ Vector3r r(s.r,s.r,s.r); mn=mn.cwise().min(s.c-r); mx=mx.cwise().max(s.c+r);}
}
Vector3r midPt() const {Vector3r mn,mx; aabb(mn,mx); return .5*(mn+mx);}
Real relDensity() const {
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp 2010-05-13 20:19:38 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp 2010-06-08 22:25:00 +0000
@@ -114,7 +114,7 @@
/* Moment Rotation Law */
Quaternionr delta( b1->state->ori * phys->initialOrientation1.conjugate() *phys->initialOrientation2 * b2->state->ori.conjugate()); //relative orientation
- AngleAxisr aa(angleAxisFromQuat(delta)); // axis of rotation - this is the Moment direction UNIT vector; angle represents the power of resistant ELASTIC moment
+ AngleAxisr aa(delta); // axis of rotation - this is the Moment direction UNIT vector; angle represents the power of resistant ELASTIC moment
if(aa.angle() > Mathr::PI) aa.angle() -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
phys->cumulativeRotation = aa.angle();
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-05-26 08:18:36 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-06-08 22:25:00 +0000
@@ -21,7 +21,7 @@
void out ( Quaternionr q )
{
- AngleAxisr aa(angleAxisFromQuat(q));
+ AngleAxisr aa(q);
std::cout << " axis: " << aa.axis()[0] << " " << aa.axis()[1] << " " << aa.axis()[2] << ", angle: " << aa.angle() << " | ";
}
@@ -123,7 +123,7 @@
delta = delta * currentContactPhysics->twistCreep;
}
- AngleAxisr aa(angleAxisFromQuat(delta)); // axis of rotation - this is the Moment direction UNIT vector; // angle represents the power of resistant ELASTIC moment
+ AngleAxisr aa(delta); // axis of rotation - this is the Moment direction UNIT vector; // angle represents the power of resistant ELASTIC moment
if (aa.angle() > Mathr::PI) aa.angle() -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
Real angle_twist(aa.angle() * aa.axis().dot(currentContactGeometry->normal));
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-06-03 03:29:25 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-06-08 22:25:00 +0000
@@ -97,7 +97,7 @@
/* Moment Rotation Law */
// NOTE this part could probably be computed in ScGeom to avoid copy/paste multiplication !!!
Quaternionr delta( b1->state->ori * phys->initialOrientation1.conjugate() *phys->initialOrientation2 * b2->state->ori.conjugate()); delta.normalize(); //relative orientation
- AngleAxisr aa(angleAxisFromQuat(delta)); // axis of rotation - this is the Moment direction UNIT vector; angle represents the power of resistant ELASTIC moment
+ AngleAxisr aa(delta); // axis of rotation - this is the Moment direction UNIT vector; angle represents the power of resistant ELASTIC moment
if(aa.angle() > Mathr::PI) aa.angle() -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
phys->cumulativeRotation = aa.angle();
=== modified file 'pkg/dem/Engine/GlobalEngine/FacetTopologyAnalyzer.cpp'
--- pkg/dem/Engine/GlobalEngine/FacetTopologyAnalyzer.cpp 2010-05-03 12:17:44 +0000
+++ pkg/dem/Engine/GlobalEngine/FacetTopologyAnalyzer.cpp 2010-06-08 22:25:00 +0000
@@ -127,7 +127,7 @@
Vector3r n1g=b1->state->ori*f1->nf, n2g=b2->state->ori*f2->nf;
//TRVAR2(n1g,n2g);
Vector3r contEdge1g=b1->state->ori*(f1->vertices[(ei+1)%3]-f1->vertices[ei]); // vector of the edge of contact in global coords
- Quaternionr q12; q12.setFromTwoVectors(n1g,(invNormals?-1.:1.)*n2g); AngleAxisr aa12(angleAxisFromQuat(q12)); Real halfAngle=.5*aa12.angle();
+ Quaternionr q12; q12.setFromTwoVectors(n1g,(invNormals?-1.:1.)*n2g); AngleAxisr aa12(q12); Real halfAngle=.5*aa12.angle();
assert(halfAngle>=0 && halfAngle<=Mathr::HALF_PI);
if(aa12.axis().dot(contEdge1g)<0 /* convex contact from the side of +n1 */ ) halfAngle*=-1.;
f1->edgeAdjHalfAngle[ei]=halfAngle;
=== modified file 'pkg/dem/Engine/GlobalEngine/MicroMacroAnalyser.cpp'
--- pkg/dem/Engine/GlobalEngine/MicroMacroAnalyser.cpp 2010-06-08 21:22:07 +0000
+++ pkg/dem/Engine/GlobalEngine/MicroMacroAnalyser.cpp 2010-06-08 22:25:00 +0000
@@ -130,7 +130,7 @@
TS.grains[Idg].sphere = CGT::Sphere(CGT::Point(pos[0],pos[1],pos[2]),rad);
// TS.grains[Idg].translation = trans;
- AngleAxisr aa(angleAxisFromQuat((*bi)->state->ori));
+ AngleAxisr aa((*bi)->state->ori);
Vector3r rotVec=aa.axis()*aa.angle();
TS.grains[Idg].rotation = CGT::Vecteur(rotVec[0],rotVec[1],rotVec[2]);
TS.box.base = CGT::Point(min(TS.box.base.x(), pos.x()-rad),
=== modified file 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp'
--- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp 2010-06-08 21:22:07 +0000
+++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp 2010-06-08 22:25:00 +0000
@@ -41,7 +41,7 @@
void NewtonIntegrator::handleClumpMemberAngAccel(Scene* scene, const body_id_t& memberId, State* memberState, State* clumpState){
// angular acceleration from: normal torque + torque generated by the force WRT particle centroid on the clump centroid
const Vector3r& m=scene->forces.getTorque(memberId); const Vector3r& f=scene->forces.getForce(memberId);
- Vector3r diffClumpAngularAccel=diagDiv(m,clumpState->inertia)+diagDiv((memberState->pos-clumpState->pos).cross(f),clumpState->inertia);
+ Vector3r diffClumpAngularAccel=m.cwise()/clumpState->inertia+(memberState->pos-clumpState->pos).cross(f).cwise()/clumpState->inertia;
// damp increment of accels on the clump, using velocities of the clump MEMBER
cundallDamp(scene->dt,m,memberState->angVel,diffClumpAngularAccel);
clumpState->angAccel+=diffClumpAngularAccel;
@@ -114,7 +114,7 @@
// rotate equation: exactAsphericalRot is disabled or the body is spherical
if (!exactAsphericalRot || (state->inertia[0]==state->inertia[1] && state->inertia[1]==state->inertia[2])){
const Vector3r& m=scene->forces.getTorque(id);
- state->angAccel=diagDiv(m,state->inertia);
+ state->angAccel=m.cwise()/state->inertia;
cundallDamp(dt,m,state->angVel,state->angAccel);
leapfrogSphericalRotate(scene,state,id,dt);
} else { // exactAsphericalRot enabled & aspherical body
@@ -144,7 +144,7 @@
leapfrogTranslate(scene,state,id,dt);
leapfrogAsphericalRotate(scene,state,id,dt,M);
} else { // exactAsphericalRot disabled or clump is spherical
- Vector3r dAngAccel=diagDiv(M,state->inertia);
+ Vector3r dAngAccel=M.cwise()/state->inertia;
cundallDamp(dt,M,state->angVel,dAngAccel);
state->angAccel+=dAngAccel;
FOREACH(Clump::memberMap::value_type mm, static_cast<Clump*>(b.get())->members){
@@ -221,12 +221,12 @@
Matrix3r A=state->ori.conjugate().toRotationMatrix(); // rotation matrix from global to local r.f.
const Vector3r l_n = state->angMom + dt/2 * M; // global angular momentum at time n
const Vector3r l_b_n = A*l_n; // local angular momentum at time n
- const Vector3r angVel_b_n = diagDiv(l_b_n,state->inertia); // local angular velocity at time n
+ const Vector3r angVel_b_n = l_b_n.cwise()/state->inertia; // local angular velocity at time n
const Quaternionr dotQ_n=DotQ(angVel_b_n,state->ori); // dQ/dt at time n
const Quaternionr Q_half = state->ori + dt/2 * dotQ_n; // Q at time n+1/2
state->angMom+=dt*M; // global angular momentum at time n+1/2
const Vector3r l_b_half = A*state->angMom; // local angular momentum at time n+1/2
- Vector3r angVel_b_half = diagDiv(l_b_half,state->inertia); // local angular velocity at time n+1/2
+ Vector3r angVel_b_half = l_b_half.cwise()/state->inertia; // local angular velocity at time n+1/2
blockRotateDOFs( state->blockedDOFs, angVel_b_half );
const Quaternionr dotQ_half=DotQ(angVel_b_half,Q_half); // dQ/dt at time n+1/2
state->ori=state->ori+dt*dotQ_half; // Q at time n+1
=== modified file 'pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp 2010-06-04 15:26:30 +0000
+++ pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp 2010-06-08 22:25:00 +0000
@@ -174,7 +174,7 @@
Quaternionr delta( de1->se3.orientation * currentContactPhysics->initialOrientation1.conjugate() *
currentContactPhysics->initialOrientation2 * de2->se3.orientation.conjugate());
- AngleAxisr aa(angleAxisFromQuat(delta)); // aa.axis() of rotation - this is the Moment direction UNIT vector; angle represents the power of resistant ELASTIC moment
+ AngleAxisr aa(delta); // aa.axis() of rotation - this is the Moment direction UNIT vector; angle represents the power of resistant ELASTIC moment
if(angle > Mathr::PI) angle -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
//This indentation is a rewrite of original equations (the two commented lines), should work exactly the same.
=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2010-06-07 17:09:08 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp 2010-06-08 22:25:00 +0000
@@ -124,11 +124,14 @@
//branch vector, FIXME : the first definition generalizes to non-spherical bodies but needs wrapped coords.
// Vector3r branch=(Body::byId(I->getId1())->state->pos-Body::byId(I->getId2())->state->pos);
Vector3r branch= gsc->normal* ( gsc->refR1+gsc->refR2 );
- // tensorial product f*branch (hand-write the tensor product to prevent matrix instanciation inside the loop by makeTensorProduct)
- stressTensor(0,0)+=f[0]*branch[0]; stressTensor(1,0)+=f[1]*branch[0]; stressTensor(2,0)+=f[2]*branch[0];
- stressTensor(0,1)+=f[0]*branch[1]; stressTensor(1,1)+=f[1]*branch[1]; stressTensor(2,1)+=f[2]*branch[1];
- stressTensor(0,2)+=f[0]*branch[2]; stressTensor(1,2)+=f[1]*branch[2]; stressTensor(2,2)+=f[2]*branch[2];
- //stressTensor+=makeTensorProduct( f, branch );
+ #if 0
+ // remove this block later
+ // tensorial product f*branch (hand-write the tensor product to prevent matrix instanciation inside the loop by makeTensorProduct)
+ //stressTensor(0,0)+=f[0]*branch[0]; stressTensor(1,0)+=f[1]*branch[0]; stressTensor(2,0)+=f[2]*branch[0];
+ //stressTensor(0,1)+=f[0]*branch[1]; stressTensor(1,1)+=f[1]*branch[1]; stressTensor(2,1)+=f[2]*branch[1];
+ //stressTensor(0,2)+=f[0]*branch[2]; stressTensor(1,2)+=f[1]*branch[2]; stressTensor(2,2)+=f[2]*branch[2];
+ #endif
+ stressTensor+=f*branch.transpose();
if( !dynCell )
{
for ( int i=0; i<3; i++ ) sumStiff[i]+=abs ( gsc->normal[i] ) *nsi->kn+ ( 1-abs ( gsc->normal[i] ) ) *nsi->ks;
=== modified file 'pkg/dem/Engine/PartialEngine/Disp2DPropLoadEngine.cpp'
--- pkg/dem/Engine/PartialEngine/Disp2DPropLoadEngine.cpp 2010-06-08 21:22:07 +0000
+++ pkg/dem/Engine/PartialEngine/Disp2DPropLoadEngine.cpp 2010-06-08 22:25:00 +0000
@@ -148,7 +148,7 @@
{
cout << "WARNING !!! your lateral boxes have not the same orientation, you're not in the case of a box imagined for creating these engines" << endl;
}
- AngleAxisr aa(angleAxisFromQuat(orientationLeftBox));
+ AngleAxisr aa(orientationLeftBox);
alpha=Mathr::PI/2.0-aa.angle(); // right if the initial orientation of the body (on the beginning of the simulation) is q =(1,0,0,0) = FromAxisAngle((0,0,1),0)
}
=== modified file 'pkg/dem/Engine/PartialEngine/KinemCNDEngine.cpp'
--- pkg/dem/Engine/PartialEngine/KinemCNDEngine.cpp 2010-05-13 20:19:38 +0000
+++ pkg/dem/Engine/PartialEngine/KinemCNDEngine.cpp 2010-06-08 22:25:00 +0000
@@ -104,7 +104,7 @@
{
cout << "WARNING !!! your lateral boxes have not the same orientation, you're not in the case of a box imagined for creating these engines" << endl;
}
- AngleAxisr aa(angleAxisFromQuat(orientationLeftBox));
+ AngleAxisr aa(orientationLeftBox);
alpha=Mathr::PI/2.0-aa.angle(); // right if the initial orientation of the body (on the beginning of the simulation) is q =(1,0,0,0) = FromAxisAngle((0,0,1),0)
}
=== modified file 'pkg/dem/Engine/PartialEngine/KinemCNLEngine.cpp'
--- pkg/dem/Engine/PartialEngine/KinemCNLEngine.cpp 2010-06-04 08:53:26 +0000
+++ pkg/dem/Engine/PartialEngine/KinemCNLEngine.cpp 2010-06-08 22:25:00 +0000
@@ -130,7 +130,7 @@
{
cout << "WARNING !!! your lateral boxes have not the same orientation, you're not in the case of a box imagined for creating these engines" << endl;
}
- AngleAxisr aa(angleAxisFromQuat(orientationLeftBox));
+ AngleAxisr aa(orientationLeftBox);
alpha=Mathr::PI/2.0-aa.angle(); // right if the initial orientation of the body (on the beginning of the simulation) is q =(1,0,0,0) = FromAxisAngle((0,0,1),0)
}
=== modified file 'pkg/dem/Engine/PartialEngine/KinemCNSEngine.cpp'
--- pkg/dem/Engine/PartialEngine/KinemCNSEngine.cpp 2010-06-04 08:53:26 +0000
+++ pkg/dem/Engine/PartialEngine/KinemCNSEngine.cpp 2010-06-08 22:25:00 +0000
@@ -114,7 +114,7 @@
{
cout << "WARNING !!! your lateral boxes have not the same orientation, you're not in the case of a box imagined for creating these engines" << endl;
}
- AngleAxisr aa(angleAxisFromQuat(orientationLeftBox));
+ AngleAxisr aa(orientationLeftBox);
alpha=Mathr::PI/2.0-aa.angle(); // right if the initial orientation of the body (on the beginning of the simulation) is q =(1,0,0,0) = FromAxisAngle((0,0,1),0)
}
=== modified file 'pkg/dem/Engine/PartialEngine/KinemCTDEngine.cpp'
--- pkg/dem/Engine/PartialEngine/KinemCTDEngine.cpp 2010-05-13 20:19:38 +0000
+++ pkg/dem/Engine/PartialEngine/KinemCTDEngine.cpp 2010-06-08 22:25:00 +0000
@@ -139,7 +139,7 @@
{
cout << "WARNING !!! your lateral boxes have not the same orientation, you're not in the case of a box imagined for creating these engines" << endl;
}
- AngleAxisr aa(angleAxisFromQuat(orientationLeftBox));
+ AngleAxisr aa(orientationLeftBox);
alpha=Mathr::PI/2.0-aa.angle(); // right if the initial orientation of the body (on the beginning of the simulation) is q =(1,0,0,0) = FromAxisAngle((0,0,1),0)
}
=== modified file 'pkg/dem/meta/ConcretePM.cpp'
--- pkg/dem/meta/ConcretePM.cpp 2010-06-08 21:22:07 +0000
+++ pkg/dem/meta/ConcretePM.cpp 2010-06-08 22:25:00 +0000
@@ -279,7 +279,7 @@
glPushMatrix();
glTranslatev(midPt);
Quaternionr q; q.setFromTwoVectors(Vector3r::UnitZ(),geom->normal);
- AngleAxisr aa(angleAxisFromQuat(q));
+ AngleAxisr aa(q);
glRotatef(aa.angle()*Mathr::RAD_TO_DEG,aa.axis()[0],aa.axis()[1],aa.axis()[2]);
glBegin(GL_POLYGON);
glColor3v(lineColor);
=== modified file 'pkg/dem/meta/Shop.cpp'
--- pkg/dem/meta/Shop.cpp 2010-06-08 21:22:07 +0000
+++ pkg/dem/meta/Shop.cpp 2010-06-08 22:25:00 +0000
@@ -192,7 +192,7 @@
FOREACH(const shared_ptr<Body>& b, *rb->bodies){
if(!b->isDynamic()) continue;
// ½(mv²+ωIω)
- ret+=.5*(b->state->mass*b->state->vel.squaredNorm()+b->state->angVel.dot(diagMult(b->state->inertia,b->state->angVel)));
+ ret+=.5*(b->state->mass*b->state->vel.squaredNorm()+b->state->angVel.dot(b->state->inertia.cwise()*b->state->angVel));
}
return ret;
}
=== modified file 'py/_eudoxos.cpp'
--- py/_eudoxos.cpp 2010-06-08 21:22:07 +0000
+++ py/_eudoxos.cpp 2010-06-08 22:25:00 +0000
@@ -152,7 +152,7 @@
if(!ge || !ph) continue;
Real r,h,theta;
boost::tie(r,h,theta)=Shop::spiralProject(ge->contactPoint,dH_dTheta,axis,NaN,theta0);
- lo=componentMinVector(lo,Vector2r(r,h)); hi=componentMaxVector(hi,Vector2r(r,h));
+ lo=lo.cwise().min(Vector2r(r,h)); hi=hi.cwise().max(Vector2r(r,h));
minD0=min(minD0,ge->refLength); maxD0=max(maxD0,ge->refLength);
minTheta=min(minTheta,theta); maxTheta=max(maxTheta,theta);
intrs.push_back(FlatInteraction(r,h,theta,i));
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2010-05-11 10:25:49 +0000
+++ py/_utils.cpp 2010-06-08 22:25:00 +0000
@@ -36,8 +36,8 @@
FOREACH(const shared_ptr<Body>& b, *Omega::instance().getScene()->bodies){
shared_ptr<Sphere> s=dynamic_pointer_cast<Sphere>(b->shape); if(!s) continue;
Vector3r rrr(s->radius,s->radius,s->radius);
- minimum=componentMinVector(minimum,Vector3r(b->state->pos-(centers?Vector3r::Zero():rrr)));
- maximum=componentMaxVector(maximum,Vector3r(b->state->pos+(centers?Vector3r::Zero():rrr)));
+ minimum=minimum.cwise().min(b->state->pos-(centers?Vector3r::Zero():rrr));
+ maximum=maximum.cwise().max(b->state->pos+(centers?Vector3r::Zero():rrr));
}
Vector3r dim=maximum-minimum;
return py::make_tuple(Vector3r(minimum+.5*cutoff*dim),Vector3r(maximum-.5*cutoff*dim));
=== modified file 'py/mathWrap/miniEigen.cpp'
--- py/mathWrap/miniEigen.cpp 2010-06-07 17:09:08 +0000
+++ py/mathWrap/miniEigen.cpp 2010-06-08 22:25:00 +0000
@@ -38,7 +38,7 @@
std::string Vector3i_str(const Vector3i & self){ return std::string("Vector3i(")+boost::lexical_cast<std::string>(self[0])+","+boost::lexical_cast<std::string>(self[1])+","+boost::lexical_cast<std::string>(self[2])+")";}
std::string Vector2r_str(const Vector2r & self){ return std::string("Vector2(")+boost::lexical_cast<std::string>(self[0])+","+boost::lexical_cast<std::string>(self[1])+")";}
std::string Vector2i_str(const Vector2i & self){ return std::string("Vector2i(")+boost::lexical_cast<std::string>(self[0])+","+boost::lexical_cast<std::string>(self[1])+")";}
-std::string Quaternionr_str(const Quaternionr & self){ AngleAxisr aa(angleAxisFromQuat(self)); return std::string("Quaternion((")+boost::lexical_cast<std::string>(aa.axis()[0])+","+boost::lexical_cast<std::string>(aa.axis()[1])+","+boost::lexical_cast<std::string>(aa.axis()[2])+"),"+boost::lexical_cast<std::string>(aa.angle())+")";}
+std::string Quaternionr_str(const Quaternionr & self){ AngleAxisr aa(self); return std::string("Quaternion((")+boost::lexical_cast<std::string>(aa.axis()[0])+","+boost::lexical_cast<std::string>(aa.axis()[1])+","+boost::lexical_cast<std::string>(aa.axis()[2])+"),"+boost::lexical_cast<std::string>(aa.angle())+")";}
std::string Matrix3r_str(const Matrix3r & self){ std::ostringstream oss; oss<<"Matrix3("; for(int i=0; i<3; i++) for(int j=0; j<3; j++) oss<<self(i,j)<<((i==2 && j==2)?")":",")<<((i<2 && j==2)?" ":""); return oss.str(); }
int Vector3r_len(){return 3;}
=== modified file 'py/pack/_packObb.cpp'
--- py/pack/_packObb.cpp 2010-05-04 13:56:05 +0000
+++ py/pack/_packObb.cpp 2010-06-08 22:25:00 +0000
@@ -24,7 +24,7 @@
Vector3r mn(inf,inf,inf), mx(-inf,-inf,-inf);
FOREACH(const Vector3r& pt, pts){
Vector3r ptT=rot*pt;
- mn=componentMinVector(mn,ptT); mx=componentMaxVector(mx,ptT);
+ mn=mn.cwise().min(ptT); mx=mx.cwise().max(ptT);
}
halfSize=.5*(mx-mn);
center=.5*(mn+mx);
=== modified file 'py/pack/_packPredicates.cpp'
--- py/pack/_packPredicates.cpp 2010-06-02 23:31:35 +0000
+++ py/pack/_packPredicates.cpp 2010-06-08 22:25:00 +0000
@@ -83,7 +83,7 @@
public:
PredicateUnion(const python::object _A, const python::object _B): PredicateBoolean(_A,_B){}
virtual bool operator()(const Vector3r& pt,Real pad) const {return obj2pred(A)(pt,pad)||obj2pred(B)(pt,pad);}
- virtual python::tuple aabb() const { Vector3r minA,maxA,minB,maxB; ttuple2vvec(obj2pred(A).aabb(),minA,maxA); ttuple2vvec(obj2pred(B).aabb(),minB,maxB); return vvec2ttuple(componentMinVector(minA,minB),componentMaxVector(maxA,maxB));}
+ virtual python::tuple aabb() const { Vector3r minA,maxA,minB,maxB; ttuple2vvec(obj2pred(A).aabb(),minA,maxA); ttuple2vvec(obj2pred(B).aabb(),minB,maxB); return vvec2ttuple(minA.cwise().min(minB),maxA.cwise().max(maxB));}
};
PredicateUnion makeUnion(const python::object& A, const python::object& B){ return PredicateUnion(A,B);}
@@ -91,7 +91,7 @@
public:
PredicateIntersection(const python::object _A, const python::object _B): PredicateBoolean(_A,_B){}
virtual bool operator()(const Vector3r& pt,Real pad) const {return obj2pred(A)(pt,pad) && obj2pred(B)(pt,pad);}
- virtual python::tuple aabb() const { Vector3r minA,maxA,minB,maxB; ttuple2vvec(obj2pred(A).aabb(),minA,maxA); ttuple2vvec(obj2pred(B).aabb(),minB,maxB); return vvec2ttuple(componentMaxVector(minA,minB),componentMinVector(maxA,maxB));}
+ virtual python::tuple aabb() const { Vector3r minA,maxA,minB,maxB; ttuple2vvec(obj2pred(A).aabb(),minA,maxA); ttuple2vvec(obj2pred(B).aabb(),minB,maxB); return vvec2ttuple(minA.cwise().max(minB),maxA.cwise().min(maxB));}
};
PredicateIntersection makeIntersection(const python::object& A, const python::object& B){ return PredicateIntersection(A,B);}
@@ -107,7 +107,7 @@
public:
PredicateSymmetricDifference(const python::object _A, const python::object _B): PredicateBoolean(_A,_B){}
virtual bool operator()(const Vector3r& pt,Real pad) const {bool inA=obj2pred(A)(pt,pad), inB=obj2pred(B)(pt,pad); return (inA && !inB) || (!inA && inB);}
- virtual python::tuple aabb() const { Vector3r minA,maxA,minB,maxB; ttuple2vvec(obj2pred(A).aabb(),minA,maxA); ttuple2vvec(obj2pred(B).aabb(),minB,maxB); return vvec2ttuple(componentMinVector(minA,minB),componentMaxVector(maxA,maxB));}
+ virtual python::tuple aabb() const { Vector3r minA,maxA,minB,maxB; ttuple2vvec(obj2pred(A).aabb(),minA,maxA); ttuple2vvec(obj2pred(B).aabb(),minB,maxB); return vvec2ttuple(minA.cwise().min(minB),maxA.cwise().max(maxB));}
};
PredicateSymmetricDifference makeSymmetricDifference(const python::object& A, const python::object& B){ return PredicateSymmetricDifference(A,B);}
@@ -158,7 +158,7 @@
sqrt((pow(A[1]-B[1],2)+pow(A[2]-B[2],2)))/ht,
sqrt((pow(A[0]-B[0],2)+pow(A[2]-B[2],2)))/ht,
sqrt((pow(A[0]-B[0],2)+pow(A[1]-B[1],2)))/ht);
- Vector3r mn=componentMinVector(A,B), mx=componentMaxVector(A,B);
+ Vector3r mn=A.cwise().min(B), mx=A.cwise().max(B);
return vvec2ttuple(mn-radius*k,mx+radius*k);
}
};
@@ -274,8 +274,8 @@
{
GtsPoint *_p=GTS_POINT(vertex);
Vector3r p(_p->x,_p->y,_p->z);
- bb->first=componentMinVector(bb->first,p);
- bb->second=componentMaxVector(bb->second,p);
+ bb->first=bb->first.cwise().min(p);
+ bb->second=bb->second.cwise().max(p);
}
/*
