yade-dev team mailing list archive

Thread
Date
[svn] r1849 - in trunk: core pkg/dem/Engine/EngineUnit

To: yade-dev@xxxxxxxxxxxxxxxxxxx
From: richefeu at BerliOS <richefeu@xxxxxxxxxxxxxxx>
Date: Wed, 8 Jul 2009 18:03:12 +0200
Author: richefeu
Date: 2009-07-08 18:03:12 +0200 (Wed, 08 Jul 2009)
New Revision: 1849

Modified:
   trunk/core/GroupRelationData.hpp
   trunk/pkg/dem/Engine/EngineUnit/BasicViscoelasticRelationships.cpp
Log:
Just try to use GroupRelationData with BasicViscolelasticRelationships.


Modified: trunk/core/GroupRelationData.hpp
===================================================================
--- trunk/core/GroupRelationData.hpp	2009-07-08 09:24:58 UTC (rev 1848)
+++ trunk/core/GroupRelationData.hpp	2009-07-08 16:03:12 UTC (rev 1849)
@@ -66,6 +66,7 @@
   ~GroupRelationData();
   
   //void activate();
+  bool isActivated() {return isActivated_;}
 
   //! Return true if the parameter exist
   //! \param name  Parameter name 

Modified: trunk/pkg/dem/Engine/EngineUnit/BasicViscoelasticRelationships.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/BasicViscoelasticRelationships.cpp	2009-07-08 09:24:58 UTC (rev 1848)
+++ trunk/pkg/dem/Engine/EngineUnit/BasicViscoelasticRelationships.cpp	2009-07-08 16:03:12 UTC (rev 1849)
@@ -12,6 +12,7 @@
 #include<yade/pkg-dem/SpheresContactGeometry.hpp>
 #include<yade/core/Omega.hpp>
 #include<yade/core/MetaBody.hpp>
+#include<yade/core/GroupRelationData.hpp>
 
 
 BasicViscoelasticRelationships::BasicViscoelasticRelationships()
@@ -30,25 +31,41 @@
 {
     if(interaction->interactionPhysics) return;
 
-    SimpleViscoelasticBodyParameters* sdec1 = static_cast<SimpleViscoelasticBodyParameters*>(b1.get());
-    SimpleViscoelasticBodyParameters* sdec2 = static_cast<SimpleViscoelasticBodyParameters*>(b2.get());
+    shared_ptr<GroupRelationData> data = (Omega::instance().getRootBody().get())->grpRelationData;
 
     interaction->interactionPhysics = shared_ptr<ViscoelasticInteraction>(new ViscoelasticInteraction());
     ViscoelasticInteraction* contactPhysics = YADE_CAST<ViscoelasticInteraction*>(interaction->interactionPhysics.get());
 
-    // Check that cn is dimensionless
-    if (sdec1->cn >= 1.0 || sdec2->cn >= 1.0) std::cout << "Warning: non dimensionless value for cn" << std::endl;
+    shared_ptr<Body> bdy1 = (*((Omega::instance().getRootBody().get())->bodies))[ interaction->getId1() ];
+    shared_ptr<Body> bdy2 = (*((Omega::instance().getRootBody().get())->bodies))[ interaction->getId2() ];
 
-    // FIXME Arbitrary ponderations... (Should be replaced by parameters that depends on groupMasks)
-    contactPhysics->kn = 0.5 * (sdec1->kn + sdec2->kn);
-    contactPhysics->ks = 0.5 * (sdec1->ks + sdec2->ks);
-    contactPhysics->cn = 0.5 * (sdec1->cn + sdec2->cn);
-    contactPhysics->cs = 0.0; // not viscosity in the tangent direction
+    if (data->isActivated())
+    {
+        int msk1 = bdy1->groupMask;
+        int msk2 = bdy2->groupMask;
+        contactPhysics->kn = (data->exists("kn")) ? data->getParameter("kn",msk1,msk2) : 0.0;
+        contactPhysics->ks = (data->exists("ks")) ? data->getParameter("ks",msk1,msk2) : 0.0;
+        contactPhysics->cn = (data->exists("cn")) ? data->getParameter("cn",msk1,msk2) : 0.0;
+        contactPhysics->cs = (data->exists("cs")) ? data->getParameter("cs",msk1,msk2) : 0.0;
+        contactPhysics->tangensOfFrictionAngle = (data->exists("mu")) ? data->getParameter("mu",msk1,msk2) : 0.0;
+    }
+    else
+    {
+        SimpleViscoelasticBodyParameters* sdec1 = static_cast<SimpleViscoelasticBodyParameters*>(b1.get());
+        SimpleViscoelasticBodyParameters* sdec2 = static_cast<SimpleViscoelasticBodyParameters*>(b2.get());
 
+        // Check that cn is dimensionless
+        if (sdec1->cn >= 1.0 || sdec2->cn >= 1.0) std::cout << "Warning: non dimensionless value for cn" << std::endl;
+    
+        // Arbitrary ponderations
+        contactPhysics->kn = 0.5 * (sdec1->kn + sdec2->kn);
+        contactPhysics->ks = 0.5 * (sdec1->ks + sdec2->ks);
+        contactPhysics->cn = 0.5 * (sdec1->cn + sdec2->cn);
+        contactPhysics->cs = 0.0; // not viscosity in the tangent direction
+        contactPhysics->tangensOfFrictionAngle = std::tan(std::min(sdec1->frictionAngle, sdec2->frictionAngle)); 
+    }
+
     // in fact, at this stage cn is a ponderation of the critical value cn_crit = 2sqrt(kn*meff)
-
-    shared_ptr<Body> bdy1 = (*((Omega::instance().getRootBody().get())->bodies))[ interaction->getId1() ];
-    shared_ptr<Body> bdy2 = (*((Omega::instance().getRootBody().get())->bodies))[ interaction->getId2() ];
     double m1,m2;
     
     if (bdy1->isClumpMember())
@@ -60,7 +77,8 @@
     }
     else
     {  
-        m1 = sdec1->mass ; 
+        ParticleParameters* pp = YADE_CAST<ParticleParameters*> ( bdy1->physicalParameters.get() );
+        m1 = pp->mass ; 
         if (!bdy1->isDynamic) m1 *= 1.0e6;
     }
     
@@ -73,17 +91,15 @@
     }
     else
     { 
-        m2 = sdec2->mass; 
+        ParticleParameters* pp = YADE_CAST<ParticleParameters*> ( bdy2->physicalParameters.get() );
+        m2 = pp->mass; 
         if (!bdy2->isDynamic) m2 *= 1.0e6;
     }
 
     contactPhysics->cn *= 2.0 * sqrt(contactPhysics->kn * ((m1 * m2) / (m1 + m2)));
  
-    contactPhysics->tangensOfFrictionAngle = std::tan(std::min(sdec1->frictionAngle, sdec2->frictionAngle)); 
-
     contactPhysics->shearForce = Vector3r(0,0,0);
     contactPhysics->prevNormal = Vector3r(0,0,0);
-
 }
 
 YADE_PLUGIN();