yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-07-29 22:36:40 +0200 (Wed, 29 Jul 2009)
New Revision: 1905

Modified:
   trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp
   trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
   trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.cpp
Log:
1. Add GenericSpheresContactClass to unit radius1, radius2 and normal between SpheresContactGeometry and Dem3DofGeom (to make it work with GlobalSitffnessTimeStepper). Partially solves https://bugs.launchpad.net/yade/+bug/399963



Modified: trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp	2009-07-29 20:05:42 UTC (rev 1904)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp	2009-07-29 20:36:40 UTC (rev 1905)
@@ -2,15 +2,26 @@
 #pragma once
 #include<yade/core/InteractionGeometry.hpp>
 
+/*! Abstract class that unites SpheresContactGeometry and Dem3DofGeom,
+	created for the purposes of GlobalStiffnessTimeStepper.
+	It exists purely on the c++, i.e. registers no attributes (the derived classes do register
+	their attributes that are initialized as references here) and doesn't exist in the yade
+	hierarchy. */
+class GenericSpheresContact: public InteractionGeometry{
+	public:
+		Vector3r normal;
+		Real refR1, refR2;
+};
+
 /*! 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)
  */
-class Dem3DofGeom: public InteractionGeometry{
+class Dem3DofGeom: public GenericSpheresContact{
 	public:
 		//! some length used to convert displacements to strains
 		Real refLength;
 		//! some unit reference vector (in the direction of the interaction)
-		Vector3r normal;
+		Vector3r &normal;
 		//! some reference point for the interaction
 		Vector3r contactPoint;
 		//! make strain go to -∞ for length going to zero
@@ -18,8 +29,10 @@
 		//! se3 of both bodies (needed to compute torque from the contact, strains etc). Must be updated at every step.
 		Se3r se31, se32;
 		//! reference radii of particles (for initial contact stiffness computation)
-		Real refR1, refR2;
+		Real &refR1, &refR2;
 
+		Dem3DofGeom(): normal(GenericSpheresContact::normal), refR1(GenericSpheresContact::refR1), refR2(GenericSpheresContact::refR2) {}
+
 		// API that needs to be implemented in derived classes
 		virtual Real displacementN();
 		virtual Vector3r displacementT(){throw;}

Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp	2009-07-29 20:05:42 UTC (rev 1904)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp	2009-07-29 20:36:40 UTC (rev 1905)
@@ -6,6 +6,7 @@
 #include<yade/core/InteractionGeometry.hpp>
 #include<yade/lib-base/yadeWm3.hpp>
 #include<yade/pkg-common/RigidBodyParameters.hpp>
+#include<yade/pkg-dem/DemXDofGeom.hpp>
 /*! Class representing geometry of two spheres in contact.
  *
  * The code under SCG_SHEAR is experimental and is used only if ElasticContactLaw::useShear is explicitly true
@@ -13,12 +14,12 @@
 
 #define SCG_SHEAR
 
-class SpheresContactGeometry: public InteractionGeometry {
+class SpheresContactGeometry: public GenericSpheresContact {
 	public:
-		Vector3r normal, // unit vector in the direction from sphere1 center to sphere2 center
-			contactPoint;
+		Vector3r& normal; // unit vector in the direction from sphere1 center to sphere2 center
+		Vector3r contactPoint;
 		Real penetrationDepth;
-		Real radius1, radius2;
+		Real &radius1, &radius2;
 
 		#ifdef SCG_SHEAR
 			//! Total value of the current shear. Update the value using updateShear.
@@ -39,7 +40,7 @@
 		#endif
 
 
-		SpheresContactGeometry():contactPoint(Vector3r::ZERO),radius1(0.),radius2(0.)
+		SpheresContactGeometry():contactPoint(Vector3r::ZERO),normal(GenericSpheresContact::normal),radius1(GenericSpheresContact::refR1),radius2(GenericSpheresContact::refR2)
 		#ifdef SCG_SHEAR
 			,shear(Vector3r::ZERO), prevNormal(Vector3r::ZERO) /*initialized to proper value by geom functor*/
 		#endif

Modified: trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.cpp	2009-07-29 20:05:42 UTC (rev 1904)
+++ trunk/pkg/dem/Engine/StandAloneEngine/GlobalStiffnessTimeStepper.cpp	2009-07-29 20:36:40 UTC (rev 1905)
@@ -9,6 +9,7 @@
 #include"GlobalStiffnessTimeStepper.hpp"
 #include<yade/pkg-dem/ElasticContactInteraction.hpp>
 #include<yade/pkg-dem/SpheresContactGeometry.hpp>
+#include<yade/pkg-dem/DemXDofGeom.hpp>
 //#include<yade/pkg-dem/MacroMicroElasticRelationships.hpp>
 #include<yade/core/Interaction.hpp>
 #include<yade/core/MetaBody.hpp>
@@ -183,10 +184,10 @@
 	FOREACH(const shared_ptr<Interaction>& contact, *rb->interactions){
 		if(!contact->isReal()) continue;
 
-		SpheresContactGeometry* geom=YADE_CAST<SpheresContactGeometry*>(contact->interactionGeometry.get()); assert(geom);
+		GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(contact->interactionGeometry.get()); assert(geom);
 		NormalShearInteraction* phys=YADE_CAST<NormalShearInteraction*>(contact->interactionPhysics.get()); assert(phys);
 		// all we need for getting stiffness
-		Vector3r& normal=geom->normal; Real& kn=phys->kn; Real& ks=phys->ks; Real& radius1=geom->radius1; Real& radius2=geom->radius2;
+		Vector3r& normal=geom->normal; Real& kn=phys->kn; Real& ks=phys->ks; Real& radius1=geom->refR1; Real& radius2=geom->refR2;
 		// FIXME? NormalShearInteraction knows nothing about whether the contact is "active" (force!=0) or not;
 		// former code: if(force==0) continue; /* disregard this interaction, it is not active */.
 		// It seems though that in such case either the interaction is accidentally at perfect equilibrium (unlikely)