yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2075
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
branch nick: trunk
timestamp: Wed 2010-03-10 21:30:46 +0100
message:
  1. NewtonIntegrator - fixes of normalizing null-vectors
  2. penetrationDepth() has been deleted from DemXDofGeom
modified:
  pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp
  pkg/dem/Engine/GlobalEngine/NewtonIntegrator.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 'pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp	2010-03-10 10:23:35 +0000
+++ pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp	2010-03-10 20:30:46 +0000
@@ -40,14 +40,6 @@
 		}
 		Vector3r strainT(){return displacementT()/refLength;}
 		Real slipToStrainTMax(Real strainTMax){return slipToDisplacementTMax(strainTMax*refLength)/refLength;}
-		
-		Real penetrationDepth(){
-			if (refR1>0) {																	//If we have both radiuses: Sphere-Sphere
-			return (refR1+refR2)-refLength-displacementN();
-			} else {																				//If we have only 1 radius: Sphere-Facet
-				return refR2-refLength-displacementN();
-			}
-		}
 
 		YADE_CLASS_BASE_DOC_ATTRS_CTOR(Dem3DofGeom,GenericSpheresContact,
 			"Abstract base class for representing contact geometry of 2 elements that has 3 degrees of freedom: normal (1 component) and shear (Vector3r, but in plane perpendicular to the normal).",

=== modified file 'pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp'
--- pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp	2010-02-14 21:23:14 +0000
+++ pkg/dem/Engine/GlobalEngine/NewtonIntegrator.cpp	2010-03-10 20:30:46 +0000
@@ -176,16 +176,30 @@
 	}
 
 }
+
 inline void NewtonIntegrator::leapfrogSphericalRotate(Scene* scene, State* state, const body_id_t& id, const Real& dt )
 {
 	blockRotateDOFs(state->blockedDOFs, state->angAccel);
 	state->angVel+=dt*state->angAccel;
-	Vector3r axis = state->angVel; Real angle = axis.Normalize();
-	Quaternionr q; q.FromAxisAngle ( axis,angle*dt );
-	state->ori = q*state->ori;
-	if(scene->forces.getMoveRotUsed() && scene->forces.getRot(id)!=Vector3r::ZERO){ Vector3r r(scene->forces.getRot(id)); Real norm=r.Normalize(); Quaternionr q; q.FromAxisAngle(r,norm); state->ori=q*state->ori; }
+	Vector3r axis = state->angVel;
+	
+	if (axis!=Vector3r::ZERO) {							//If we have an angular velocity, we make a rotation
+		Quaternionr q;
+		Real angle = axis.Normalize();
+		q.FromAxisAngle ( axis,angle*dt );
+		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.Normalize();
+		q.FromAxisAngle(r,norm);
+		state->ori=q*state->ori;
+	}
 	state->ori.Normalize();
 }
+
 void NewtonIntegrator::leapfrogAsphericalRotate(Scene* scene, State* state, const body_id_t& id, const Real& dt, const Vector3r& M){
 	Matrix3r A; state->ori.Conjugate().ToRotationMatrix(A); // rotation matrix from global to local r.f.
 	const Vector3r l_n = state->angMom + dt/2 * M; // global angular momentum at time n
@@ -200,7 +214,14 @@
 	const Quaternionr dotQ_half=DotQ(angVel_b_half,Q_half); // dQ/dt at time n+1/2
 	state->ori+=dt*dotQ_half; // Q at time n+1
 	state->angVel=state->ori.Rotate(angVel_b_half); // global angular velocity at time n+1/2
-	if(scene->forces.getMoveRotUsed() && scene->forces.getRot(id)!=Vector3r::ZERO){ Vector3r r(scene->forces.getRot(id)); Real norm=r.Normalize(); Quaternionr q; q.FromAxisAngle(r,norm); state->ori=q*state->ori; }
+
+	if(scene->forces.getMoveRotUsed() && scene->forces.getRot(id)!=Vector3r::ZERO) {
+		Vector3r r(scene->forces.getRot(id));
+		Real norm=r.Normalize();
+		Quaternionr q;
+		q.FromAxisAngle(r,norm);
+		state->ori=q*state->ori;
+	}
 	state->ori.Normalize(); 
 }