yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2106
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Fri 2010-03-26 21:39:30 +0100
message:
  - Some documentation and code cleaning.
modified:
  pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp
  pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.hpp


--
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/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp'
--- pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp	2010-01-10 09:09:32 +0000
+++ pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp	2010-03-26 20:39:30 +0000
@@ -16,68 +16,40 @@
 
 
 
-void Ip2_FrictMat_FrictMat_FrictPhys::go(	  const shared_ptr<Material>& b1
+void Ip2_FrictMat_FrictMat_FrictPhys::go( const shared_ptr<Material>& b1
 					, const shared_ptr<Material>& b2
 					, const shared_ptr<Interaction>& interaction)
 {
-	
-	//ScGeom* interactionGeometry = YADE_CAST<ScGeom*>(interaction->interactionGeometry.get());
-	
-	//if(interactionGeometry)
+	if(!interaction->interactionPhysics)
 	{
-		if(!interaction->interactionPhysics)
-		{
-			const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1);
-			const shared_ptr<FrictMat>& mat2 = YADE_PTR_CAST<FrictMat>(b2);
-			
-			if (!interaction->interactionPhysics) interaction->interactionPhysics = shared_ptr<FrictPhys>(new FrictPhys());
-			
-			const shared_ptr<FrictPhys>& contactPhysics = YADE_PTR_CAST<FrictPhys>(interaction->interactionPhysics);
-
-			Real Ea 	= mat1->young;
-			Real Eb 	= mat2->young;
-			Real Va 	= mat1->poisson;
-			Real Vb 	= mat2->poisson;
-			#if 0
-				Real Da 	= interactionGeometry->radius1; // FIXME - multiply by factor of sphere interaction distance (so sphere interacts at bigger range that its geometrical size)
-				Real Db 	= interactionGeometry->radius2; // FIXME - as above
-				Vector3r normal=interactionGeometry->normal;
-			#else
-				Real Da,Db; Vector3r normal;
-				ScGeom* scg=dynamic_cast<ScGeom*>(interaction->interactionGeometry.get());
-				Dem3DofGeom* d3dg=dynamic_cast<Dem3DofGeom*>(interaction->interactionGeometry.get());
-				if(scg){ Da=scg->radius1; Db=scg->radius2; normal=scg->normal; }
-				else if(d3dg){Da=d3dg->refR1>0?d3dg->refR1:2*d3dg->refR2; Db=d3dg->refR2>0?d3dg->refR2:d3dg->refR1; normal=d3dg->normal; }
-				else throw runtime_error("Ip2_FrictMat_FrictMat_FrictPhys: geometry is neither ScGeom nor Dem3DofGeom");
-			#endif
-			
-			Real fa 	= mat1->frictionAngle;
-			Real fb 	= mat2->frictionAngle;
-
-			//Real Eab	= 2*Ea*Eb/(Ea+Eb);
-			//Real Vab	= 2*Va*Vb/(Va+Vb);
-
-			Real Dinit 	= Da+Db; 			// FIXME - is it just a sum?
-			//Real Sinit 	= Mathr::PI * std::pow( std::min(Da,Db) , 2);
-
-			Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
-			Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);//harmonic average of two stiffnesses with ks=V*kn for each sphere
-
-
-			contactPhysics->initialKn			= Kn;
-			contactPhysics->initialKs			= Ks;
-//cerr << "Ks: " <<       contactPhysics->initialKs			<< endl;
-			contactPhysics->frictionAngle			= std::min(fa,fb); // FIXME - this is actually a waste of memory space, just like initialKs and initialKn
-			contactPhysics->tangensOfFrictionAngle		= std::tan(contactPhysics->frictionAngle); 
-
-			contactPhysics->prevNormal 			= normal;
-			contactPhysics->initialEquilibriumDistance	= Dinit;			
-
-			contactPhysics->kn = contactPhysics->initialKn;
-			contactPhysics->ks = contactPhysics->initialKs;
-			contactPhysics->equilibriumDistance = contactPhysics->initialEquilibriumDistance;
-
-		}	
+		const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1);
+		const shared_ptr<FrictMat>& mat2 = YADE_PTR_CAST<FrictMat>(b2);
+		if (!interaction->interactionPhysics)
+			interaction->interactionPhysics = shared_ptr<FrictPhys>(new FrictPhys());
+		const shared_ptr<FrictPhys>& contactPhysics = YADE_PTR_CAST<FrictPhys>(interaction->interactionPhysics);
+
+		Real Ea 	= mat1->young;
+		Real Eb 	= mat2->young;
+		Real Va 	= mat1->poisson;
+		Real Vb 	= mat2->poisson;
+		
+		Real Da,Db; Vector3r normal;
+		//FIXME : dynamic casts here???!!!
+		ScGeom* scg=dynamic_cast<ScGeom*>(interaction->interactionGeometry.get());
+		Dem3DofGeom* d3dg=dynamic_cast<Dem3DofGeom*>(interaction->interactionGeometry.get());
+		if(scg){Da=scg->radius1; Db=scg->radius2; normal=scg->normal;}
+		else if(d3dg){Da=d3dg->refR1>0?d3dg->refR1:2*d3dg->refR2; Db=d3dg->refR2>0?d3dg->refR2:d3dg->refR1; normal=d3dg->normal;}
+		else throw runtime_error("Ip2_FrictMat_FrictMat_FrictPhys: geometry is neither ScGeom nor Dem3DofGeom");
+		//harmonic average of the two stiffnesses when (Di.Ei/2) is the stiffness of sphere "i"
+		Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);
+		//same for shear stiffness
+		Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);
+		
+		contactPhysics->frictionAngle = std::min(mat1->frictionAngle,mat2->frictionAngle);
+		contactPhysics->tangensOfFrictionAngle = std::tan(contactPhysics->frictionAngle); 
+		contactPhysics->prevNormal = normal;
+		contactPhysics->kn = Kn;
+		contactPhysics->ks = Ks;
 		return;
 	}
 	throw runtime_error("Ip2_FrictMat_FrictMat_FrictPhys currently fails for non-ScGeom geometry!");

=== modified file 'pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.hpp'
--- pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.hpp	2010-02-18 21:08:54 +0000
+++ pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.hpp	2010-03-26 20:39:30 +0000
@@ -16,7 +16,7 @@
 			const shared_ptr<Material>& b2,
 			const shared_ptr<Interaction>& interaction);
 	FUNCTOR2D(FrictMat,FrictMat);
-	YADE_CLASS_BASE_DOC(Ip2_FrictMat_FrictMat_FrictPhys,InteractionPhysicsFunctor,"Create a :yref:`FrictPhys` from two :yref:`FrictMats<FrictMat>`. Most parameters are averaged, but the exact algorithm is not documented. Only interactions with :yref:`ScGeom` or :yref:`Dem3DofGeom` geometry are meaningfully accepted; run-time typecheck can make this functor unnecessarily slow in general. Such design is problematic in itself, though -- from http://www.mail-archive.com/yade-dev@xxxxxxxxxxxxxxxxxxx/msg02603.html:\n\n\t\t\tYou have to suppose some exact type of InteractionGeometry in the Ip2 functor, but you don't know anything about it (Ip2 only guarantees you get certain InteractionPhysics types, via the dispatch mechanism).\n\n\t\t\tThat means, unless you use Ig2 functor producing the desired type, the code will break (crash or whatever). The right behavior would be either to accept any type (what we have now, at least in principle), or really enforce InteractionGeometry type of the interation passed to that particular Ip2 functor.\n\nEtc.");
+	YADE_CLASS_BASE_DOC(Ip2_FrictMat_FrictMat_FrictPhys,InteractionPhysicsFunctor,"Create a :yref:`FrictPhys` from two :yref:`FrictMats<FrictMat>`. The compliance of one sphere under symetric point loads is defined in here as 1/(E.D), with E the stifness of the sphere and D its diameter, and corresponds to a compliance 1/(2.E.D) from each side. The compliance of the contact itself will be the sum of compliances from each sphere, i.e. 1/(2.E.D1)+1/(2.E.D1) in the general case, or 1/(E.D) in the special case of equal sizes. Note that summing compliances corresponds to an harmonic average of stiffnesss, which is how kn is actually computed in the :yref:Ip2_FrictMat_FrictMat_FrictPhys functor. The friction angle of the contact is defined as the minimum angle of the two materials in contact. \n\n\ Only interactions with :yref:`ScGeom` or :yref:`Dem3DofGeom` geometry are meaningfully accepted; run-time typecheck can make this functor unnecessarily slow in general. Such design is problematic in itself, though -- from http://www.mail-archive.com/yade-dev@xxxxxxxxxxxxxxxxxxx/msg02603.html:\n\n\t\t\tYou have to suppose some exact type of InteractionGeometry in the Ip2 functor, but you don't know anything about it (Ip2 only guarantees you get certain InteractionPhysics types, via the dispatch mechanism).\n\n\t\t\tThat means, unless you use Ig2 functor producing the desired type, the code will break (crash or whatever). The right behavior would be either to accept any type (what we have now, at least in principle), or really enforce InteractionGeometry type of the interation passed to that particular Ip2 functor.\n\nEtc.");
 };
 REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_FrictPhys);