本文整理汇总了C++中Vector3r::dot方法的典型用法代码示例。如果您正苦于以下问题：C++ Vector3r::dot方法的具体用法？C++ Vector3r::dot怎么用？C++ Vector3r::dot使用的例子？那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类Vector3r的用法示例。
在下文中一共展示了Vector3r::dot方法的13个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: operator
bool operator()(const Vector3r& pt, Real pad=0.) const {
Real u=(pt.dot(c12)-c1.dot(c12))/(ht*ht); // normalized coordinate along the c1--c2 axis
if((u*ht<0+pad) || (u*ht>ht-pad)) return false; // out of cylinder along the axis
Real axisDist=((pt-c1).cross(pt-c2)).norm()/ht;
if(axisDist>radius-pad) return false;
return true;
}示例2:
notInNotch(const Vector3r& _c, const Vector3r& _edge, const Vector3r& _normal, Real _aperture){
c=_c;
edge=_edge; edge.normalize();
normal=_normal; normal-=edge*edge.dot(normal); normal.normalize();
inside=edge.cross(normal);
aperture=_aperture;
// LOG_DEBUG("edge="<<edge<<", normal="<<normal<<", inside="<<inside<<", aperture="<<aperture);
}示例3: a
bool Law2_ScGeom_ElectrostaticPhys::go(shared_ptr<IGeom>& iGeom, shared_ptr<IPhys>& iPhys, Interaction* interaction)
{
// Physic
ElectrostaticPhys* phys=static_cast<ElectrostaticPhys*>(iPhys.get());
// Geometry
ScGeom* geom=static_cast<ScGeom*>(iGeom.get());
// Get bodies properties
Body::id_t id1 = interaction->getId1();
Body::id_t id2 = interaction->getId2();
const shared_ptr<Body> b1 = Body::byId(id1,scene);
const shared_ptr<Body> b2 = Body::byId(id2,scene);
State* s1 = b1->state.get();
State* s2 = b2->state.get();
// geometric parameters
Real a((geom->radius1+geom->radius2)/2.);
bool isNew=false;
// Speeds
Vector3r shiftVel=scene->isPeriodic ? Vector3r(scene->cell->velGrad*scene->cell->hSize*interaction->cellDist.cast<Real>()) : Vector3r::Zero();
Vector3r shift2 = scene->isPeriodic ? Vector3r(scene->cell->hSize*interaction->cellDist.cast<Real>()): Vector3r::Zero();
Vector3r relV = geom->getIncidentVel(s1, s2, scene->dt, shift2, shiftVel, false );
Real undot = relV.dot(geom->normal); // Normal velocity norm
if(-geom->penetrationDepth > a && -geom->penetrationDepth > 10.*phys->DebyeLength) { return false; }
// inititalization
if(phys->u == -1. ) {phys->u = -geom->penetrationDepth; isNew=true;}
// Solve normal
normalForce_DL_Adim(phys,geom, undot,isNew);
// Solve shear and torques
Vector3r C1 = Vector3r::Zero();
Vector3r C2 = Vector3r::Zero();
computeShearForceAndTorques_log(phys, geom, s1, s2, C1, C2);
// Apply!
scene->forces.addForce(id1,phys->normalForce+phys->shearForce);
scene->forces.addTorque(id1,C1);
scene->forces.addForce(id2,-(phys->normalForce+phys->shearForce));
scene->forces.addTorque(id2,C2);
return true;
}示例4: ScGridCoGeom
//! O/
bool Ig2_Sphere_GridConnection_ScGridCoGeom::go( const shared_ptr<Shape>& cm1,
const shared_ptr<Shape>& cm2,
const State& state1, const State& state2, const Vector3r& shift2, const bool& force,
const shared_ptr<Interaction>& c)
{ // Useful variables :
const State* sphereSt = YADE_CAST<const State*>(&state1);
//const State* gridCoSt = YADE_CAST<const State*>(&state2);
Sphere* sphere = YADE_CAST<Sphere*>(cm1.get());
GridConnection* gridCo = YADE_CAST<GridConnection*>(cm2.get());
//GridNode* gridNo1 = YADE_CAST<GridNode*>(gridCo->node1->shape.get());
//GridNode* gridNo2 = YADE_CAST<GridNode*>(gridCo->node2->shape.get());
State* gridNo1St = YADE_CAST<State*>(gridCo->node1->state.get());
State* gridNo2St = YADE_CAST<State*>(gridCo->node2->state.get());
bool isNew = !c->geom;
shared_ptr<ScGridCoGeom> scm;
if (!isNew) scm = YADE_PTR_CAST<ScGridCoGeom>(c->geom);
Vector3r segt = gridCo->getSegment();
float len = gridCo->getLength();
Vector3r branch = sphereSt->pos - gridNo1St->pos;
float relPos = branch.dot(segt)/(len*len);
Vector3r fictiousPos=gridNo1St->pos+relPos*segt;
Vector3r branchP = fictiousPos - sphereSt->pos;
float dist = branchP.norm();
if(isNew){
if(dist > (sphere->radius + gridCo->radius)) return false;
else {scm=shared_ptr<ScGridCoGeom>(new ScGridCoGeom()); c->geom=scm;}
}
if(dist <= (sphere->radius + gridCo->radius)){
scm->refR1=sphere->radius; //FIXME don't know why I have to do that ...
scm->refR2=gridCo->radius;
scm->id3=gridCo->node1->getId();
scm->id4=gridCo->node2->getId();
scm->relPos=relPos;
Vector3r normal=branchP/dist;
scm->penetrationDepth=sphere->radius+gridCo->radius-dist;
scm->fictiousState.pos = gridNo1St->pos+segt*relPos;
scm->fictiousState.vel = (1-relPos)*gridNo1St->vel + relPos*gridNo2St->vel;
scm->fictiousState.angVel =
((1-relPos)*gridNo1St->angVel + relPos*gridNo2St->angVel).dot(segt/len)*segt/len //twist part : interpolated
+ segt.cross(gridNo2St->vel - gridNo1St->vel);// non-twist part : defined from nodes velocities
scm->contactPoint = sphereSt->pos+normal*(sphere->radius-0.5*scm->penetrationDepth);
scm->precompute(state1,scm->fictiousState,scene,c,normal,isNew,shift2,true);//use sphere-sphere precompute (with a virtual sphere)
}
return true;
}示例5: action
void TorqueRecorder::action(){
totalTorque=0;
Vector3r tmpAxis = rotationAxis.normalized();
FOREACH(Body::id_t id, ids){
if (!(scene->bodies->exists(id))) continue;
Body* b=Body::byId(id,scene).get();
Vector3r tmpPos = b->state->pos;
Vector3r radiusVector = tmpAxis.cross(tmpAxis.cross(tmpPos - zeroPoint));
totalTorque+=tmpAxis.dot(scene->forces.getTorque(id)+radiusVector.cross(scene->forces.getForce(id)));
};
//Save data to a file
out<<scene->iter<<" "<<totalTorque<<"\n";
out.close();
}示例6: if
bool Law2_ScGeom_MindlinPhys_Mindlin::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
const Real& dt = scene->dt; // get time step
Body::id_t id1 = contact->getId1(); // get id body 1
Body::id_t id2 = contact->getId2(); // get id body 2
State* de1 = Body::byId(id1,scene)->state.get();
State* de2 = Body::byId(id2,scene)->state.get();
ScGeom* scg = static_cast<ScGeom*>(ig.get());
MindlinPhys* phys = static_cast<MindlinPhys*>(ip.get());
const shared_ptr<Body>& b1=Body::byId(id1,scene);
const shared_ptr<Body>& b2=Body::byId(id2,scene);
bool useDamping=(phys->betan!=0. || phys->betas!=0. || phys->alpha!=0.);
bool LinDamp=true;
if (phys->alpha!=0.) {LinDamp=false;} // use non linear damping
// tangential and normal stiffness coefficients, recomputed from betan,betas at every step
Real cn=0, cs=0;
/****************/
/* NORMAL FORCE */
/****************/
Real uN = scg->penetrationDepth; // get overlapping
if (uN<0) {
if (neverErase) {phys->shearForce = phys->normalForce = Vector3r::Zero(); phys->kn=phys->ks=0; return true;}
else return false;
}
/* Hertz-Mindlin's formulation (PFC)
Note that the normal stiffness here is a secant value (so as it is cannot be used in the GSTS)
In the first place we get the normal force and then we store kn to be passed to the GSTS */
Real Fn = phys->kno*std::pow(uN,1.5); // normal Force (scalar)
if (includeAdhesion) {
Fn -= phys->adhesionForce; // include adhesion force to account for the effect of Van der Waals interactions
phys->isAdhesive = (Fn<0); // set true the bool to count the number of adhesive contacts
}
phys->normalForce = Fn*scg->normal; // normal Force (vector)
if (calcEnergy){
Real R=scg->radius1*scg->radius2/(scg->radius1+scg->radius2);
phys->radius=pow((Fn+(includeAdhesion?phys->adhesionForce:0.))*pow(R,3/2.)/phys->kno,1/3.); // attribute not used anywhere, we do not need it
}
/*******************************/
/* TANGENTIAL NORMAL STIFFNESS */
/*******************************/
phys->kn = 3./2.*phys->kno*std::pow(uN,0.5); // here we store the value of kn to compute the time step
/******************************/
/* TANGENTIAL SHEAR STIFFNESS */
/******************************/
phys->ks = phys->kso*std::pow(uN,0.5); // get tangential stiffness (this is a tangent value, so we can pass it to the GSTS)
/************************/
/* DAMPING COEFFICIENTS */
/************************/
// Inclusion of local damping if requested
// viscous damping is defined for both linear and non-linear elastic case
if (useDamping && LinDamp){
Real mbar = (!b1->isDynamic() && b2->isDynamic()) ? de2->mass : ((!b2->isDynamic() && b1->isDynamic()) ? de1->mass : (de1->mass*de2->mass / (de1->mass + de2->mass))); // get equivalent mass if both bodies are dynamic, if not set it equal to the one of the dynamic body
//Real mbar = de1->mass*de2->mass / (de1->mass + de2->mass); // equivalent mass
Real Cn_crit = 2.*sqrt(mbar*phys->kn); // Critical damping coefficient (normal direction)
Real Cs_crit = 2.*sqrt(mbar*phys->ks); // Critical damping coefficient (shear direction)
// Note: to compare with the analytical solution you provide cn and cs directly (since here we used a different method to define c_crit)
cn = Cn_crit*phys->betan; // Damping normal coefficient
cs = Cs_crit*phys->betas; // Damping tangential coefficient
if(phys->kn<0 || phys->ks<0){ cerr<<"Negative stiffness kn="<<phys->kn<<" ks="<<phys->ks<<" for ##"<<b1->getId()<<"+"<<b2->getId()<<", step "<<scene->iter<<endl; }
}
else if (useDamping){ // (see Tsuji, 1992)
Real mbar = (!b1->isDynamic() && b2->isDynamic()) ? de2->mass : ((!b2->isDynamic() && b1->isDynamic()) ? de1->mass : (de1->mass*de2->mass / (de1->mass + de2->mass))); // get equivalent mass if both bodies are dynamic, if not set it equal to the one of the dynamic body
cn = phys->alpha*sqrt(mbar)*pow(uN,0.25); // normal viscous coefficient, see also [Antypov2011] eq. 10
cs = cn; // same value for shear viscous coefficient
}
/***************/
/* SHEAR FORCE */
/***************/
Vector3r& shearElastic = phys->shearElastic; // reference for shearElastic force
// Define shifts to handle periodicity
const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero();
const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(contact->cellDist): Vector3r::Zero();
// 1. Rotate shear force
shearElastic = scg->rotate(shearElastic);
Vector3r prev_FsElastic = shearElastic; // save shear force at previous time step
// 2. Get incident velocity, get shear and normal components
Vector3r incidentV = scg->getIncidentVel(de1, de2, dt, shift2, shiftVel, preventGranularRatcheting);
Vector3r incidentVn = scg->normal.dot(incidentV)*scg->normal; // contact normal velocity
Vector3r incidentVs = incidentV - incidentVn; // contact shear velocity
// 3. Get shear force (incrementally)
shearElastic = shearElastic - phys->ks*(incidentVs*dt);
/**************************************/
/* VISCOUS DAMPING (Normal direction) */
//.........这里部分代码省略.........
示例7:
/*! @brief Recalculates inertia tensor of a body after translation away from (default) or towards its centroid.
*
* @param I inertia tensor in the original coordinates; it is assumed to be upper-triangular (elements below the diagonal are ignored).
* @param m mass of the body; if positive, translation is away from the centroid; if negative, towards centroid.
* @param off offset of the new origin from the original origin
* @return inertia tensor in the new coordinate system; the matrix is symmetric.
*/
Matrix3r woo::Volumetric::inertiaTensorTranslate(const Matrix3r& I, const Real m, const Vector3r& off){
// short eigen implementation; check it gives the same result as above
return I+m*(off.dot(off)*Matrix3r::Identity()-off*off.transpose());
}示例8: if
bool Ig2_Facet_Sphere_ScGeom::go(const shared_ptr<Shape>& cm1,
const shared_ptr<Shape>& cm2,
const State& state1,
const State& state2,
const Vector3r& shift2,
const bool& force,
const shared_ptr<Interaction>& c)
{
TIMING_DELTAS_START();
const Se3r& se31=state1.se3; const Se3r& se32=state2.se3;
Facet* facet = static_cast<Facet*>(cm1.get());
/* could be written as (needs to be tested):
* Vector3r cl=se31.orientation.Conjugate()*(se32.position-se31.position);
*/
Matrix3r facetAxisT=se31.orientation.toRotationMatrix();
Matrix3r facetAxis = facetAxisT.transpose();
// local orientation
Vector3r cl = facetAxis*(se32.position + shift2 - se31.position); // "contact line" in facet-local coords
//
// BEGIN everything in facet-local coordinates
//
Vector3r normal = facet->normal;
Real L = normal.dot(cl);
if (L<0) {normal=-normal; L=-L; }
Real sphereRadius = static_cast<Sphere*>(cm2.get())->radius;
if (L>sphereRadius && !c->isReal() && !force) { // no contact, but only if there was no previous contact; ortherwise, the constitutive law is responsible for setting Interaction::isReal=false
TIMING_DELTAS_CHECKPOINT("Ig2_Facet_Sphere_ScGeom");
return false;
}
Vector3r cp = cl - L*normal;
const Vector3r* ne = facet->ne;
Real penetrationDepth=0;
Real bm = ne[0].dot(cp);
int m=0;
for (int i=1; i<3; ++i)
{
Real b=ne[i].dot(cp);
if (bm<b) {bm=b; m=i;}
}
Real sh = sphereRadius*shrinkFactor;
Real icr = facet->icr-sh;
if (icr<0)
{
LOG_WARN("a radius of a facet's inscribed circle less than zero! So, shrinkFactor is too large and would be reduced to zero.");
shrinkFactor=0;
icr = facet->icr;
sh = 0;
}
if (bm<icr) // contact with facet's surface
{
penetrationDepth = sphereRadius - L;
normal.normalize();
}
else
{
cp = cp + ne[m]*(icr-bm);
if (cp.dot(ne[(m-1<0)?2:m-1])>icr) // contact with vertex m
// cp = facet->vertices[m];
cp = facet->vu[m]*(facet->vl[m]-sh);
else if (cp.dot(ne[m=(m+1>2)?0:m+1])>icr) // contact with vertex m+1
// cp = facet->vertices[(m+1>2)?0:m+1];
cp = facet->vu[m]*(facet->vl[m]-sh);
normal = cl-cp;
Real norm=normal.norm(); normal/=norm;
penetrationDepth = sphereRadius - norm;
}
//
// END everything in facet-local coordinates
//
if (penetrationDepth>0 || c->isReal())
{
shared_ptr<ScGeom> scm;
bool isNew = !c->geom;
if (c->geom)
scm = YADE_PTR_CAST<ScGeom>(c->geom);
else
scm = shared_ptr<ScGeom>(new ScGeom());
normal = facetAxisT*normal; // in global orientation
scm->contactPoint = se32.position + shift2 - (sphereRadius-0.5*penetrationDepth)*normal;
scm->penetrationDepth = penetrationDepth;
scm->radius1 = 2*sphereRadius;
scm->radius2 = sphereRadius;
if (isNew) c->geom = scm;
scm->precompute(state1,state2,scene,c,normal,isNew,shift2,false/*avoidGranularRatcheting only for sphere-sphere*/);
TIMING_DELTAS_CHECKPOINT("Ig2_Facet_Sphere_ScGeom");
return true;
}
//.........这里部分代码省略.........
示例9: inertiaTranslate
Matrix3r inertiaTranslate(const Matrix3r& I, const Real m, const Vector3r& off){
return I+m*(off.dot(off)*Matrix3r::Identity()-off*off.transpose());
}示例10: action
void RockLiningGlobal::action(){
const double PI = 3.14159;
if (openingCreated == true && installed == false){
double angleInterval = 2.0*PI/static_cast<double>(totalNodes);
for (int n=0; n<totalNodes;n++){
double currentAngle = 0.0 + n*angleInterval; /* from 0 degrees east */
double unitX = cos(currentAngle);
double unitY = sin(currentAngle);
Vector3r searchDir(unitX,0,unitY);
vector<double> distanceFrOpening; vector<int> IDs;
double outerRadius = openingRad + 1.0;
FOREACH(const shared_ptr<Body>& b, *scene->bodies){
if (!b) continue;
if (b->isClump() == true) continue;
PotentialBlock* pb=static_cast<PotentialBlock*>(b->shape.get());
if(!pb) continue;
if(pb->isBoundary == true || pb->erase== true || pb->isLining==true){continue;}
State* state1 = b->state.get();
Vector3r intersectionPt(0,0,0);
if ( installLining(pb, state1, startingPoint, searchDir, outerRadius, intersectionPt )){
IDs.push_back(b->id);
distanceFrOpening.push_back((intersectionPt-startingPoint).norm());
//std::cout<<"currentAngle: "<<currentAngle<<", b->id: "<<b->id<<", dist: "<<(intersectionPt-startingPoint).norm()<<endl;
}
}
/* find closest block */
int totalBlocks = IDs.size();
double closestDistance = 100000.0;
int closestID=0;
for (int i=0; i<totalBlocks; i++){
if ( distanceFrOpening[i] < closestDistance){
closestID = IDs[i];
closestDistance = distanceFrOpening[i];
}
}
stickIDs.push_back(closestID);
IDs.clear();distanceFrOpening.clear();
//std::cout<<"closestID: "<<closestID<<endl;
/* find intersection with edges of polygon */
Vector3r jointIntersection (0,0,0);
State* state1 = Body::byId(closestID,scene)->state.get();
Shape* shape1 = Body::byId(closestID,scene)->shape.get();
PotentialBlock *pb=static_cast<PotentialBlock*>(shape1);
int totalPlanes = pb->a.size();
int intersectNo = 0;
Vector3r nodeLocalPos(0,0,0); Vector3r nodeGlobalPos(0,0,0);
//std::cout<<"totalPlanes: "<<totalPlanes<<endl;
double closestPlaneDist = 1000000;
for (int i=0; i<totalPlanes; i++){
Vector3r plane = state1->ori*Vector3r(pb->a[i], pb->b[i], pb->c[i]); double planeD = plane.dot(state1->pos) + pb->d[i] +pb->r;
if ( intersectPlane(pb, state1,startingPoint,searchDir, outerRadius, jointIntersection, plane, planeD)){
double distance = jointIntersection.norm();
if (distance < closestPlaneDist){
closestPlaneDist = distance;
nodeLocalPos = state1->ori.conjugate()*(jointIntersection-state1->pos);
nodeGlobalPos=jointIntersection;
}
}
}
if(nodeGlobalPos.norm() > 1.03*openingRad){ nodeGlobalPos=1.03*openingRad*searchDir;}
//if(nodeGlobalPos.norm() < 0.98*openingRad){ continue;}
//initOverlap = interfaceTension/interfaceStiffness;
nodeGlobalPos = nodeGlobalPos + searchDir*initOverlap;
localCoordinates.push_back(nodeLocalPos);
refPos.push_back(nodeGlobalPos);
int nodeID = insertNode(nodeGlobalPos, lumpedMass, contactLength);
blockIDs.push_back(nodeID); //(nodeID); //(closestID);
refOri.push_back(Quaternionr::Identity()); //(state1->ori);
installed = true;
axialForces.push_back(0.0);
shearForces.push_back(0.0);
moment.push_back(0.0);
sigmaMax.push_back(0.0);
sigmaMin.push_back(0.0);
displacement.push_back(0.0);
radialDisplacement.push_back(0.0);
}
totalNodes=blockIDs.size();
/* Assembling global stiffness matrix */
for (int n=0; n<totalNodes;n++){
int nextID = n+1;
if(nextID==totalNodes){nextID=0;}
double Length = (refPos[nextID]-refPos[n]).norm();
lengthNode.push_back(Length);
Vector3r localDir = refPos[nextID]-refPos[n];
localDir.normalize(); refDir.push_back(localDir);
double angle = acos(localDir.dot(Vector3r(1,0,0)));
Vector3r signAngle = Vector3r(1,0.0,0).cross(localDir);
if (signAngle.dot(Vector3r(0,-1.0,0)) < 0.0){angle =2.0*PI - angle;}
refAngle.push_back(angle);
//.........这里部分代码省略.........
示例11: recompute
void SphereClumpGeom::recompute(int _div, bool failOk, bool fastOnly){
if((centers.empty() && radii.empty()) || centers.size()!=radii.size()){
if(failOk) { makeInvalid(); return;}
throw std::runtime_error("SphereClumpGeom.recompute: centers and radii must have the same length (len(centers)="+to_string(centers.size())+", len(radii)="+to_string(radii.size())+"), and may not be empty.");
}
div=_div;
// one single sphere: simple
if(centers.size()==1){
pos=centers[0];
ori=Quaternionr::Identity();
volume=(4/3.)*M_PI*pow(radii[0],3);
inertia=Vector3r::Constant((2/5.)*volume*pow(radii[0],2));
equivRad=radii[0];
return;
}
volume=0;
Vector3r Sg=Vector3r::Zero();
Matrix3r Ig=Matrix3r::Zero();
if(_div<=0){
// non-intersecting: Steiner's theorem
for(size_t i=0; i<centers.size(); i++){
const Real& r(radii[i]); const Vector3r& x(centers[i]);
Real v=(4/3.)*M_PI*pow(r,3);
volume+=v;
Sg+=v*x;
Ig+=woo::Volumetric::inertiaTensorTranslate(Vector3r::Constant((2/5.)*v*pow(r,2)).asDiagonal(),v,-1.*x);
}
} else {
// intersecting: grid sampling
Real rMin=Inf; AlignedBox3r aabb;
for(size_t i=0; i<centers.size(); i++){
aabb.extend(centers[i]+Vector3r::Constant(radii[i]));
aabb.extend(centers[i]-Vector3r::Constant(radii[i]));
rMin=min(rMin,radii[i]);
}
if(rMin<=0){
if(failOk){ makeInvalid(); return; }
throw std::runtime_error("SphereClumpGeom.recompute: minimum radius must be positive (not "+to_string(rMin)+")");
}
Real dx=rMin/_div; Real dv=pow(dx,3);
long nCellsApprox=(aabb.sizes()/dx).prod();
// don't compute anything, it would take too long
if(fastOnly && nCellsApprox>1e5){ makeInvalid(); return; }
if(nCellsApprox>1e8) LOG_WARN("SphereClumpGeom: space grid has "<<nCellsApprox<<" cells, computing inertia can take a long time.");
Vector3r x;
for(x.x()=aabb.min().x()+dx/2.; x.x()<aabb.max().x(); x.x()+=dx){
for(x.y()=aabb.min().y()+dx/2.; x.y()<aabb.max().y(); x.y()+=dx){
for(x.z()=aabb.min().z()+dx/2.; x.z()<aabb.max().z(); x.z()+=dx){
for(size_t i=0; i<centers.size(); i++){
if((x-centers[i]).squaredNorm()<pow(radii[i],2)){
volume+=dv;
Sg+=dv*x;
Ig+=dv*(x.dot(x)*Matrix3r::Identity()-x*x.transpose())+/*along princial axes of dv; perhaps negligible?*/Matrix3r(Vector3r::Constant(dv*pow(dx,2)/6.).asDiagonal());
break;
}
}
}
}
}
}
woo::Volumetric::computePrincipalAxes(volume,Sg,Ig,pos,ori,inertia);
equivRad=(inertia.array()/volume).sqrt().mean(); // mean of radii of gyration
}示例12:
void Law2_ScGeom_CapillaryPhys_Capillarity::checkFusion()
{
//Reset fusion numbers
InteractionContainer::iterator ii = scene->interactions->begin();
InteractionContainer::iterator iiEnd = scene->interactions->end();
for( ; ii!=iiEnd ; ++ii ) {
if ((*ii)->isReal()) {
if (!hertzOn) static_cast<CapillaryPhys*>((*ii)->phys.get())->fusionNumber=0;
else static_cast<MindlinCapillaryPhys*>((*ii)->phys.get())->fusionNumber=0;
}
}
list< shared_ptr<Interaction> >::iterator firstMeniscus, lastMeniscus, currentMeniscus;
Real angle1 = -1.0; Real angle2 = -1.0;
for ( int i=0; i< bodiesMenisciiList.size(); ++i ) { // i is the index (or id) of the body being tested
CapillaryPhys* cundallInteractionPhysics1=NULL;
MindlinCapillaryPhys* mindlinInteractionPhysics1=NULL;
CapillaryPhys* cundallInteractionPhysics2=NULL;
MindlinCapillaryPhys* mindlinInteractionPhysics2=NULL;
if ( !bodiesMenisciiList[i].empty() ) {
lastMeniscus = bodiesMenisciiList[i].end();
for ( firstMeniscus=bodiesMenisciiList[i].begin(); firstMeniscus!=lastMeniscus; ++firstMeniscus ) { //FOR EACH MENISCUS ON THIS BODY...
currentMeniscus = firstMeniscus; ++currentMeniscus;
if (!hertzOn) {
cundallInteractionPhysics1 = YADE_CAST<CapillaryPhys*>((*firstMeniscus)->phys.get());
if (i == (*firstMeniscus)->getId1()) angle1=cundallInteractionPhysics1->Delta1;//get angle of meniscus1 on body i
else angle1=cundallInteractionPhysics1->Delta2;
}
else {
mindlinInteractionPhysics1 = YADE_CAST<MindlinCapillaryPhys*>((*firstMeniscus)->phys.get());
if (i == (*firstMeniscus)->getId1()) angle1=mindlinInteractionPhysics1->Delta1;//get angle of meniscus1 on body i
else angle1=mindlinInteractionPhysics1->Delta2;
}
for ( ;currentMeniscus!= lastMeniscus; ++currentMeniscus) {//... CHECK FUSION WITH ALL OTHER MENISCII ON THE BODY
if (!hertzOn) {
cundallInteractionPhysics2 = YADE_CAST<CapillaryPhys*>((*currentMeniscus)->phys.get());
if (i == (*currentMeniscus)->getId1()) angle2=cundallInteractionPhysics2->Delta1;//get angle of meniscus2 on body i
else angle2=cundallInteractionPhysics2->Delta2;
}
else {
mindlinInteractionPhysics2 = YADE_CAST<MindlinCapillaryPhys*>((*currentMeniscus)->phys.get());
if (i == (*currentMeniscus)->getId1()) angle2=mindlinInteractionPhysics2->Delta1;//get angle of meniscus2 on body i
else angle2=mindlinInteractionPhysics2->Delta2;
}
if (angle1==0 || angle2==0) cerr << "THIS SHOULD NOT HAPPEN!!"<< endl;
//cerr << "angle1 = " << angle1 << " | angle2 = " << angle2 << endl;
Vector3r normalFirstMeniscus = YADE_CAST<ScGeom*>((*firstMeniscus)->geom.get())->normal;
Vector3r normalCurrentMeniscus = YADE_CAST<ScGeom*>((*currentMeniscus)->geom.get())->normal;
Real normalDot = 0;
if ((*firstMeniscus)->getId1() == (*currentMeniscus)->getId1() || (*firstMeniscus)->getId2() == (*currentMeniscus)->getId2()) normalDot = normalFirstMeniscus.dot(normalCurrentMeniscus);
else normalDot = - (normalFirstMeniscus.dot(normalCurrentMeniscus));
Real normalAngle = 0;
if (normalDot >= 0 ) normalAngle = Mathr::FastInvCos0(normalDot);
else normalAngle = ((Mathr::PI) - Mathr::FastInvCos0(-(normalDot)));
if ((angle1+angle2)*Mathr::DEG_TO_RAD > normalAngle) {
if (!hertzOn) {++(cundallInteractionPhysics1->fusionNumber); ++(cundallInteractionPhysics2->fusionNumber);}//count +1 if 2 meniscii are overlaping
else {++(mindlinInteractionPhysics1->fusionNumber); ++(mindlinInteractionPhysics2->fusionNumber);}
};
}
}
}
}
}示例13: if
//! O/
bool Ig2_Sphere_GridConnection_ScGridCoGeom::go( const shared_ptr<Shape>& cm1,
const shared_ptr<Shape>& cm2,
const State& state1, const State& state2, const Vector3r& shift2, const bool& force,
const shared_ptr<Interaction>& c)
{ // Useful variables :
const State* sphereSt = YADE_CAST<const State*>(&state1);
//const State* gridCoSt = YADE_CAST<const State*>(&state2);
Sphere* sphere = YADE_CAST<Sphere*>(cm1.get());
GridConnection* gridCo = YADE_CAST<GridConnection*>(cm2.get());
GridNode* gridNo1 = YADE_CAST<GridNode*>(gridCo->node1->shape.get());
GridNode* gridNo2 = YADE_CAST<GridNode*>(gridCo->node2->shape.get());
State* gridNo1St = YADE_CAST<State*>(gridCo->node1->state.get());
State* gridNo2St = YADE_CAST<State*>(gridCo->node2->state.get());
bool isNew = !c->geom;
shared_ptr<ScGridCoGeom> scm;
if (!isNew) scm = YADE_PTR_CAST<ScGridCoGeom>(c->geom);
else {scm = shared_ptr<ScGridCoGeom>(new ScGridCoGeom());}
Vector3r segt = gridCo->getSegment();
Real len = gridCo->getLength();
Vector3r spherePos = sphereSt->pos - shift2;
Vector3r branch = spherePos - gridNo1St->pos;
Vector3r branchN = spherePos - gridNo2St->pos;
for(int i=0;i<3;i++){
if(abs(branch[i])<1e-14) branch[i]=0.0;
if(abs(branchN[i])<1e-14) branchN[i]=0.0;
}
Real relPos = branch.dot(segt)/(len*len);
if(scm->isDuplicate==2 && scm->trueInt!=c->id2)return true; //the contact will be deleted into the Law.
scm->isDuplicate=0;
scm->trueInt=-1;
if(relPos<=0){ // if the sphere projection is BEFORE the segment ...
if(gridNo1->ConnList.size()>1){// if the node is not an extremity of the Grid (only one connection)
for(int unsigned i=0;i<gridNo1->ConnList.size();i++){ // ... loop on all the Connections of the same Node ...
GridConnection* GC = (GridConnection*)gridNo1->ConnList[i]->shape.get();
if(GC==gridCo)continue;// self comparison.
Vector3r segtCandidate1 = GC->node1->state->pos - gridNo1St->pos; // (be sure of the direction of segtPrev to compare relPosPrev.)
Vector3r segtCandidate2 = GC->node2->state->pos - gridNo1St->pos;
Vector3r segtPrev = segtCandidate1.norm()>segtCandidate2.norm() ? segtCandidate1:segtCandidate2;
for(int j=0;j<3;j++){
if(abs(segtPrev[j])<1e-14) segtPrev[j]=0.0;
}
Real relPosPrev = (branch.dot(segtPrev))/(segtPrev.norm()*segtPrev.norm());
// ... and check whether the sphere projection is before the neighbours connections too.
const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,gridNo1->ConnList[i]->getId());
if(relPosPrev<=0){ //if the sphere projection is outside both the current Connection AND this neighbouring connection, then create the interaction if the neighbour did not already do it before.
if(intr && intr->isReal() && isNew) return false;
if(intr && intr->isReal() && !isNew) {scm->isDuplicate=1;/*cout<<"Declare "<<c->id1<<"-"<<c->id2<<" as duplicated."<<endl;*/}
}
else{//the sphere projection is outside the current Connection but inside the previous neighbour. The contact has to be handled by the Prev GridConnection, not here.
if (isNew)return false;
else {
//cout<<"The contact "<<c->id1<<"-"<<c->id2<<" HAVE to be copied and deleted NOW."<<endl ;
scm->isDuplicate=1;
scm->trueInt=-1;
return true;}
}
}
}
}
//Exactly the same but in the case the sphere projection is AFTER the segment.
else if(relPos>=1){
if(gridNo2->ConnList.size()>1){
for(int unsigned i=0;i<gridNo2->ConnList.size();i++){
GridConnection* GC = (GridConnection*)gridNo2->ConnList[i]->shape.get();
if(GC==gridCo)continue;// self comparison.
Vector3r segtCandidate1 = GC->node1->state->pos - gridNo2St->pos;
Vector3r segtCandidate2 = GC->node2->state->pos - gridNo2St->pos;
Vector3r segtNext = segtCandidate1.norm()>segtCandidate2.norm() ? segtCandidate1:segtCandidate2;
for(int j=0;j<3;j++){
if(abs(segtNext[j])<1e-14) segtNext[j]=0.0;
}
Real relPosNext = (branchN.dot(segtNext))/(segtNext.norm()*segtNext.norm());
const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,gridNo2->ConnList[i]->getId());
if(relPosNext<=0){ //if the sphere projection is outside both the current Connection AND this neighbouring connection, then create the interaction if the neighbour did not already do it before.
if(intr && intr->isReal() && isNew) return false;
if(intr && intr->isReal() && !isNew) {scm->isDuplicate=1;/*cout<<"Declare "<<c->id1<<"-"<<c->id2<<" as duplicated."<<endl;*/}
}
else{//the sphere projection is outside the current Connection but inside the previous neighbour. The contact has to be handled by the Prev GridConnection, not here.
if (isNew)return false;
else {//cout<<"The contact "<<c->id1<<"-"<<c->id2<<" HAVE to be copied and deleted NOW."<<endl ;
scm->isDuplicate=1 ;
scm->trueInt=-1 ;
return true;}
}
}
}
}
else if (isNew && relPos<=0.5){
if(gridNo1->ConnList.size()>1){// if the node is not an extremity of the Grid (only one connection)
for(int unsigned i=0;i<gridNo1->ConnList.size();i++){ // ... loop on all the Connections of the same Node ...
GridConnection* GC = (GridConnection*)gridNo1->ConnList[i]->shape.get();
if(GC==gridCo)continue;// self comparison.
Vector3r segtCandidate1 = GC->node1->state->pos - gridNo1St->pos; // (be sure of the direction of segtPrev to compare relPosPrev.)
Vector3r segtCandidate2 = GC->node2->state->pos - gridNo1St->pos;
Vector3r segtPrev = segtCandidate1.norm()>segtCandidate2.norm() ? segtCandidate1:segtCandidate2;
for(int j=0;j<3;j++){
//.........这里部分代码省略.........