yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01127
[svn] r1718 - in trunk: examples examples/dynamic_simulation_tests pkg/common/DataClass/PhysicalParameters pkg/common/Engine/DeusExMachina pkg/dem pkg/dem/DataClass/InteractionGeometry pkg/dem/Engine/EngineUnit pkg/dem/Engine/StandAloneEngine pkg/dem/PreProcessor
Author: sega
Date: 2009-03-12 18:27:35 +0100 (Thu, 12 Mar 2009)
New Revision: 1718
Added:
trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp
trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp
Removed:
trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp
trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp
trunk/pkg/dem/Engine/StandAloneEngine/SimpleViscoelasticContactLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/SimpleViscoelasticContactLaw.hpp
Modified:
trunk/examples/STLImporterTest.py
trunk/examples/dynamic_simulation_tests/ringSimpleViscoelastic.py
trunk/pkg/common/DataClass/PhysicalParameters/RigidBodyParameters.cpp
trunk/pkg/common/DataClass/PhysicalParameters/RigidBodyParameters.hpp
trunk/pkg/common/Engine/DeusExMachina/RotationEngine.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
trunk/pkg/dem/PreProcessor/MembraneTest.cpp
trunk/pkg/dem/PreProcessor/TestSimpleViscoelastic.cpp
trunk/pkg/dem/SConscript
Log:
1. Remove zeroPoint
2. Add prefix ef2_ to Spheres_Viscoelastic_SimpleViscoelasticContactLaw
Modified: trunk/examples/STLImporterTest.py
===================================================================
--- trunk/examples/STLImporterTest.py 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/examples/STLImporterTest.py 2009-03-12 17:27:35 UTC (rev 1718)
@@ -64,7 +64,7 @@
## Create physical information about the interaction.
MetaEngine('InteractionPhysicsMetaEngine',[EngineUnit('SimpleViscoelasticRelationships')]),
## Constitutive law
- MetaEngine('ConstitutiveLawDispatcher',[EngineUnit('Spheres_Viscoelastic_SimpleViscoelasticContactLaw')]),
+ MetaEngine('ConstitutiveLawDispatcher',[EngineUnit('ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw')]),
## Apply gravity
DeusExMachina('GravityEngine',{'gravity':[0,-9.81,0]}),
## Cundall damping must been disabled!
Modified: trunk/examples/dynamic_simulation_tests/ringSimpleViscoelastic.py
===================================================================
--- trunk/examples/dynamic_simulation_tests/ringSimpleViscoelastic.py 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/examples/dynamic_simulation_tests/ringSimpleViscoelastic.py 2009-03-12 17:27:35 UTC (rev 1718)
@@ -68,9 +68,7 @@
MetaEngine('InteractionPhysicsMetaEngine',[EngineUnit('SimpleViscoelasticRelationships')]),
## Constitutive law
- MetaEngine('ConstitutiveLawDispatcher',[EngineUnit('Spheres_Viscoelastic_SimpleViscoelasticContactLaw')]),
- ## or InteractionSolver
- #StandAloneEngine('SimpleViscoelasticContactLaw'),
+ MetaEngine('ConstitutiveLawDispatcher',[EngineUnit('ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw')]),
DeusExMachina('GravityEngine',{'gravity':[0,-9.81,0]}),
Modified: trunk/pkg/common/DataClass/PhysicalParameters/RigidBodyParameters.cpp
===================================================================
--- trunk/pkg/common/DataClass/PhysicalParameters/RigidBodyParameters.cpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/common/DataClass/PhysicalParameters/RigidBodyParameters.cpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -16,7 +16,6 @@
acceleration = Vector3r(0,0,0);
angularAcceleration = Vector3r(0,0,0);
angularVelocity=Vector3r(0,0,0);
- zeroPoint=Vector3r(0,0,0);
}
RigidBodyParameters::~RigidBodyParameters()
@@ -28,7 +27,6 @@
ParticleParameters::registerAttributes();
REGISTER_ATTRIBUTE(inertia);
REGISTER_ATTRIBUTE(angularVelocity);
-// REGISTER_ATTRIBUTE(zeroPoint); // must be set by kinematic Engine
}
YADE_PLUGIN();
Modified: trunk/pkg/common/DataClass/PhysicalParameters/RigidBodyParameters.hpp
===================================================================
--- trunk/pkg/common/DataClass/PhysicalParameters/RigidBodyParameters.hpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/common/DataClass/PhysicalParameters/RigidBodyParameters.hpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -20,10 +20,6 @@
// state
,angularAcceleration
,angularVelocity;
-
- /// It is the rotation center of kinematic (non-isDymanic) body.
- /// It is non-serializable and must be set by kinematic Engine.
- Vector3r zeroPoint;
RigidBodyParameters ();
virtual ~RigidBodyParameters ();
Modified: trunk/pkg/common/Engine/DeusExMachina/RotationEngine.cpp
===================================================================
--- trunk/pkg/common/Engine/DeusExMachina/RotationEngine.cpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/common/Engine/DeusExMachina/RotationEngine.cpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -74,7 +74,6 @@
std::vector<int>::const_iterator iiEnd = subscribedBodies.end();
Real dt = Omega::instance().getTimeStep();
- // time = dt;
Quaternionr q;
q.FromAxisAngle(rotationAxis,angularVelocity*dt);
@@ -86,16 +85,19 @@
{
RigidBodyParameters * rb = static_cast<RigidBodyParameters*>((*bodies)[*ii]->physicalParameters.get());
+ rb->angularVelocity = rotationAxis*angularVelocity;
+
if(rotateAroundZero)
- rb->se3.position = q*(rb->se3.position-zeroPoint)+zeroPoint; // for RotatingBox
+ {
+ const Vector3r l = rb->se3.position-zeroPoint;
+ rb->se3.position = q*l+zeroPoint;
+ rb->velocity = rb->angularVelocity.Cross(l);
+ }
rb->se3.orientation = q*rb->se3.orientation;
-
rb->se3.orientation.Normalize();
rb->se3.orientation.ToAxisAngle(ax,an);
- rb->angularVelocity = rotationAxis*angularVelocity;
- rb->velocity = Vector3r(0,0,0);
}
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -9,7 +9,7 @@
// At least one virtual method must be in the .cpp file (!!!)
SpheresContactGeometry::~SpheresContactGeometry(){};
-void SpheresContactGeometry::updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const RigidBodyParameters* rbp1, const RigidBodyParameters* rbp2, bool isDynamic1, bool isDynamic2, Real dt, bool avoidGranularRatcheting){
+void SpheresContactGeometry::updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const RigidBodyParameters* rbp1, const RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting){
Vector3r axis;
Real angle;
@@ -17,13 +17,7 @@
// approximated rotations
axis = prevNormal.Cross(normal);
shearForce -= shearForce.Cross(axis);
- //angle = dt*0.5*currentContactGeometry->normal.Dot(de1->angularVelocity+de2->angularVelocity);
- //FIXME: if one body is kinematic then assumed its rotation centre does not lies along the normal
- //(i.e. virtual sphere, which replaces this kinematic body on contact, does not rotate)
- Vector3r summaryAngularVelocity(Vector3r::ZERO);
- if (isDynamic1) summaryAngularVelocity += rbp1->angularVelocity;
- if (isDynamic2) summaryAngularVelocity += rbp2->angularVelocity;
- angle = dt*0.5*normal.Dot(summaryAngularVelocity);
+ angle = dt*0.5*normal.Dot(rbp1->angularVelocity + rbp2->angularVelocity);
axis = angle*normal;
shearForce -= shearForce.Cross(axis);
@@ -48,10 +42,13 @@
* (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN,
* "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena",
* ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors) */
- c1x = (isDynamic1) ? radius1*normal : x - rbp1->zeroPoint;
- c2x = (isDynamic2) ? -radius2*normal : x - rbp2->zeroPoint;
+
+ // FIXME: For sphere-facet contact this will give an erroneous value of relative velocity...
+ c1x = radius1*normal;
+ c2x = -radius2*normal;
}
else {
+ // FIXME: It is correct for sphere-sphere and sphere-facet contact
c1x = (x - rbp1->se3.position);
c2x = (x - rbp2->se3.position);
}
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -100,7 +100,7 @@
SpheresContactGeometry():contactPoint(Vector3r::ZERO),radius1(0),radius2(0),facetContactFace(0.),hasShear(false),pos1(Vector3r::ZERO),pos2(Vector3r::ZERO),ori1(Quaternionr::IDENTITY),ori2(Quaternionr::IDENTITY),cp1rel(Quaternionr::IDENTITY),cp2rel(Quaternionr::IDENTITY),d1(0),d2(0),d0(0),initRelOri12(Quaternionr::IDENTITY){createIndex();}
virtual ~SpheresContactGeometry();
- void updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const RigidBodyParameters* rbp1, const RigidBodyParameters* rbp2, bool isDynamic1, bool isDynamic2, Real dt, bool avoidGranularRatcheting=true);
+ void updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const RigidBodyParameters* rbp1, const RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting=true);
REGISTER_ATTRIBUTES(/* no attributes from base class */,
(normal)
Deleted: trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -1,74 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2009 by Sergei Dorofeenko *
-* sega@xxxxxxxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include"Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp"
-#include<yade/pkg-dem/SpheresContactGeometry.hpp>
-#include<yade/pkg-dem/ViscoelasticInteraction.hpp>
-#include<yade/pkg-common/RigidBodyParameters.hpp>
-YADE_PLUGIN("Spheres_Viscoelastic_SimpleViscoelasticContactLaw");
-
-void Spheres_Viscoelastic_SimpleViscoelasticContactLaw::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody){
-
- SpheresContactGeometry* geom=static_cast<SpheresContactGeometry*>(_geom.get());
- ViscoelasticInteraction* phys=static_cast<ViscoelasticInteraction*>(_phys.get());
-
- int id1 = I->getId1();
- int id2 = I->getId2();
-
- shared_ptr<BodyContainer>& bodies = rootBody->bodies;
-
- RigidBodyParameters* de1 = YADE_CAST<RigidBodyParameters*>((*bodies)[id1]->physicalParameters.get());
- RigidBodyParameters* de2 = YADE_CAST<RigidBodyParameters*>((*bodies)[id2]->physicalParameters.get());
-
- bool isDynamic1 = (*bodies)[id1]->isDynamic;
- bool isDynamic2 = (*bodies)[id2]->isDynamic;
-
- Vector3r& shearForce = phys->shearForce;
- if (I->isNew) shearForce=Vector3r(0,0,0);
-
- Real dt = Omega::instance().getTimeStep();
-
- Vector3r axis = phys->prevNormal.Cross(geom->normal);
- shearForce -= shearForce.Cross(axis);
- Vector3r summaryAngularVelocity(0,0,0);
- if (isDynamic1) summaryAngularVelocity += de1->angularVelocity;
- if (isDynamic2) summaryAngularVelocity += de2->angularVelocity;
- Real angle = dt*0.5*geom->normal.Dot(summaryAngularVelocity);
- axis = angle*geom->normal;
- shearForce -= shearForce.Cross(axis);
-
- Vector3r x = geom->contactPoint;
- Vector3r c1x = (x - de1->se3.position);
- Vector3r c2x = (x - de2->se3.position);
- /// The following definition of c1x and c2x is to avoid "granular ratcheting"
- /// (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN,
- /// "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena",
- /// ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors)
- Vector3r _c1x_ = (isDynamic1) ? geom->radius1*geom->normal : x - de1->zeroPoint;
- Vector3r _c2x_ = (isDynamic2) ? -geom->radius2*geom->normal : x - de2->zeroPoint;
- Vector3r relativeVelocity = (de2->velocity+de2->angularVelocity.Cross(_c2x_)) - (de1->velocity+de1->angularVelocity.Cross(_c1x_));
- Real normalVelocity = geom->normal.Dot(relativeVelocity);
- Vector3r shearVelocity = relativeVelocity-normalVelocity*geom->normal;
- shearForce -= (phys->ks*dt+phys->cs)*shearVelocity;
-
- phys->normalForce = ( phys->kn * std::max(geom->penetrationDepth,(Real) 0) - phys->cn * normalVelocity ) * geom->normal;
- phys->prevNormal = geom->normal;
-
- Real maxFs = phys->normalForce.SquaredLength() * std::pow(phys->tangensOfFrictionAngle,2);
- if( shearForce.SquaredLength() > maxFs )
- {
- maxFs = Mathr::Sqrt(maxFs) / shearForce.Length();
- shearForce *= maxFs;
- }
-
- Vector3r f = phys->normalForce + shearForce;
- addForce (id1,-f,rootBody);
- addForce (id2, f,rootBody);
- addTorque(id1,-c1x.Cross(f),rootBody);
- addTorque(id2, c2x.Cross(f),rootBody);
-}
Deleted: trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -1,24 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2009 by Sergei Dorofeenko *
-* sega@xxxxxxxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg-common/ConstitutiveLaw.hpp>
-
-/// This class provides linear viscoelastic contact model
-class Spheres_Viscoelastic_SimpleViscoelasticContactLaw : public ConstitutiveLaw
-{
- public :
- virtual void go(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*, MetaBody*);
- NEEDS_BEX("Force","Momentum");
- FUNCTOR2D(SpheresContactGeometry,ViscoelasticInteraction);
- REGISTER_CLASS_AND_BASE(Spheres_Viscoelastic_SimpleViscoelasticContactLaw,ConstitutiveLaw);
-};
-REGISTER_SERIALIZABLE(Spheres_Viscoelastic_SimpleViscoelasticContactLaw);
-
-
Copied: trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp (from rev 1717, trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp)
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -0,0 +1,69 @@
+/*************************************************************************
+* Copyright (C) 2009 by Sergei Dorofeenko *
+* sega@xxxxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+#include"ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp"
+#include<yade/pkg-dem/SpheresContactGeometry.hpp>
+#include<yade/pkg-dem/ViscoelasticInteraction.hpp>
+#include<yade/pkg-common/RigidBodyParameters.hpp>
+YADE_PLUGIN("ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw");
+
+void ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody){
+
+ SpheresContactGeometry* geom=static_cast<SpheresContactGeometry*>(_geom.get());
+ ViscoelasticInteraction* phys=static_cast<ViscoelasticInteraction*>(_phys.get());
+
+ int id1 = I->getId1();
+ int id2 = I->getId2();
+
+ shared_ptr<BodyContainer>& bodies = rootBody->bodies;
+
+ RigidBodyParameters* de1 = YADE_CAST<RigidBodyParameters*>((*bodies)[id1]->physicalParameters.get());
+ RigidBodyParameters* de2 = YADE_CAST<RigidBodyParameters*>((*bodies)[id2]->physicalParameters.get());
+
+ Vector3r& shearForce = phys->shearForce;
+ if (I->isNew) shearForce=Vector3r(0,0,0);
+
+ Real dt = Omega::instance().getTimeStep();
+
+ Vector3r axis = phys->prevNormal.Cross(geom->normal);
+ shearForce -= shearForce.Cross(axis);
+ Real angle = dt*0.5*geom->normal.Dot(de1->angularVelocity + de2->angularVelocity);
+ axis = angle*geom->normal;
+ shearForce -= shearForce.Cross(axis);
+
+ Vector3r x = geom->contactPoint;
+ Vector3r c1x = (x - de1->se3.position);
+ Vector3r c2x = (x - de2->se3.position);
+ /// The following definition of c1x and c2x is to avoid "granular ratcheting"
+ /// (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN,
+ /// "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena",
+ /// ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors)
+ //Vector3r _c1x_ = geom->radius1*geom->normal;
+ //Vector3r _c2x_ = -geom->radius2*geom->normal;
+ //Vector3r relativeVelocity = (de2->velocity+de2->angularVelocity.Cross(_c2x_)) - (de1->velocity+de1->angularVelocity.Cross(_c1x_));
+ Vector3r relativeVelocity = (de2->velocity+de2->angularVelocity.Cross(c2x)) - (de1->velocity+de1->angularVelocity.Cross(c1x));
+ Real normalVelocity = geom->normal.Dot(relativeVelocity);
+ Vector3r shearVelocity = relativeVelocity-normalVelocity*geom->normal;
+ shearForce -= (phys->ks*dt+phys->cs)*shearVelocity;
+
+ phys->normalForce = ( phys->kn * std::max(geom->penetrationDepth,(Real) 0) - phys->cn * normalVelocity ) * geom->normal;
+ phys->prevNormal = geom->normal;
+
+ Real maxFs = phys->normalForce.SquaredLength() * std::pow(phys->tangensOfFrictionAngle,2);
+ if( shearForce.SquaredLength() > maxFs )
+ {
+ maxFs = Mathr::Sqrt(maxFs) / shearForce.Length();
+ shearForce *= maxFs;
+ }
+
+ Vector3r f = phys->normalForce + shearForce;
+ addForce (id1,-f,rootBody);
+ addForce (id2, f,rootBody);
+ addTorque(id1,-c1x.Cross(f),rootBody);
+ addTorque(id2, c2x.Cross(f),rootBody);
+}
Property changes on: trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp
___________________________________________________________________
Name: svn:mergeinfo
+
Copied: trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp (from rev 1717, trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp)
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -0,0 +1,24 @@
+/*************************************************************************
+* Copyright (C) 2009 by Sergei Dorofeenko *
+* sega@xxxxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+#pragma once
+
+#include<yade/pkg-common/ConstitutiveLaw.hpp>
+
+/// This class provides linear viscoelastic contact model
+class ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw : public ConstitutiveLaw
+{
+ public :
+ virtual void go(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*, MetaBody*);
+ NEEDS_BEX("Force","Momentum");
+ FUNCTOR2D(SpheresContactGeometry,ViscoelasticInteraction);
+ REGISTER_CLASS_AND_BASE(ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw,ConstitutiveLaw);
+};
+REGISTER_SERIALIZABLE(ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw);
+
+
Property changes on: trunk/pkg/dem/Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp
___________________________________________________________________
Name: svn:mergeinfo
+
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -104,9 +104,6 @@
BodyMacroParameters* de1 = YADE_CAST<BodyMacroParameters*>((*bodies)[id1]->physicalParameters.get());
BodyMacroParameters* de2 = YADE_CAST<BodyMacroParameters*>((*bodies)[id2]->physicalParameters.get());
- bool isDynamic1 = (*bodies)[id1]->isDynamic;
- bool isDynamic2 = (*bodies)[id2]->isDynamic;
-
Vector3r& shearForce = currentContactPhysics->shearForce;
if (contact->isNew) shearForce=Vector3r(0,0,0);
@@ -122,18 +119,9 @@
Real angle;
/// Here is the code with approximated rotations ///
-
-
-
axis = currentContactPhysics->prevNormal.Cross(currentContactGeometry->normal);
shearForce -= shearForce.Cross(axis);
- //angle = dt*0.5*currentContactGeometry->normal.Dot(de1->angularVelocity+de2->angularVelocity);
- //FIXME: if one body is kinematic then assumed its rotation centre does not lies along the normal
- //(i.e. virtual sphere, which replaces this kinematic body on contact, does not rotate)
- Vector3r summaryAngularVelocity(0,0,0);
- if (isDynamic1) summaryAngularVelocity += de1->angularVelocity;
- if (isDynamic2) summaryAngularVelocity += de2->angularVelocity;
- angle = dt*0.5*currentContactGeometry->normal.Dot(summaryAngularVelocity);
+ angle = dt*0.5*currentContactGeometry->normal.Dot(de1->angularVelocity+de2->angularVelocity);
axis = angle*currentContactGeometry->normal;
shearForce -= shearForce.Cross(axis);
@@ -161,9 +149,10 @@
/// (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN,
/// "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena",
/// ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors)
- Vector3r _c1x_ = (isDynamic1) ? currentContactGeometry->radius1*currentContactGeometry->normal : x - de1->zeroPoint;
- Vector3r _c2x_ = (isDynamic2) ? -currentContactGeometry->radius2*currentContactGeometry->normal : x - de2->zeroPoint;
- Vector3r relativeVelocity = (de2->velocity+de2->angularVelocity.Cross(_c2x_)) - (de1->velocity+de1->angularVelocity.Cross(_c1x_));
+ //Vector3r _c1x_ = currentContactGeometry->radius1*currentContactGeometry->normal;
+ //Vector3r _c2x_ = -currentContactGeometry->radius2*currentContactGeometry->normal;
+ //Vector3r relativeVelocity = (de2->velocity+de2->angularVelocity.Cross(_c2x_)) - (de1->velocity+de1->angularVelocity.Cross(_c1x_));
+ Vector3r relativeVelocity = (de2->velocity+de2->angularVelocity.Cross(c2x)) - (de1->velocity+de1->angularVelocity.Cross(c1x));
Vector3r shearVelocity = relativeVelocity-currentContactGeometry->normal.Dot(relativeVelocity)*currentContactGeometry->normal;
Vector3r shearDisplacement = shearVelocity*dt;
shearForce -= currentContactPhysics->ks*shearDisplacement;
Deleted: trunk/pkg/dem/Engine/StandAloneEngine/SimpleViscoelasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/SimpleViscoelasticContactLaw.cpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/dem/Engine/StandAloneEngine/SimpleViscoelasticContactLaw.cpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -1,153 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2004 by Olivier Galizzi *
-* olivier.galizzi@xxxxxxx *
-* *
-* Copyright (C) 2008 by Sergei Dorofeenko *
-* sega@xxxxxxxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include"SimpleViscoelasticContactLaw.hpp"
-#include<yade/pkg-dem/SimpleViscoelasticBodyParameters.hpp>
-#include<yade/pkg-dem/ViscoelasticInteraction.hpp>
-#include<yade/pkg-dem/SpheresContactGeometry.hpp>
-#include<yade/core/Omega.hpp>
-#include<yade/core/MetaBody.hpp>
-#include<yade/pkg-common/Force.hpp>
-#include<yade/pkg-common/Momentum.hpp>
-#include<yade/core/PhysicalAction.hpp>
-
-
-SimpleViscoelasticContactLaw::SimpleViscoelasticContactLaw() : InteractionSolver() , actionForce(new Force) , actionMomentum(new Momentum)
-{
-// sdecGroupMask=1;
- momentRotationLaw = true;
- actionForceIndex = actionForce->getClassIndex();
- actionMomentumIndex = actionMomentum->getClassIndex();
-}
-
-
-void SimpleViscoelasticContactLaw::registerAttributes()
-{
- InteractionSolver::registerAttributes();
-// REGISTER_ATTRIBUTE(sdecGroupMask);
- REGISTER_ATTRIBUTE(momentRotationLaw);
-}
-
-
-void SimpleViscoelasticContactLaw::action(MetaBody* ncb)
-{
- shared_ptr<BodyContainer>& bodies = ncb->bodies;
-
- Real dt = Omega::instance().getTimeStep();
-
-/// Non Permanents Links ///
-
- InteractionContainer::iterator ii = ncb->transientInteractions->begin();
- InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
- for( ; ii!=iiEnd ; ++ii )
- {
- if ((*ii)->isReal)
- {
- const shared_ptr<Interaction>& contact = *ii;
- int id1 = contact->getId1();
- int id2 = contact->getId2();
-
- //if( !( (*bodies)[id1]->getGroupMask() & (*bodies)[id2]->getGroupMask() & sdecGroupMask) ) continue;
-
- SpheresContactGeometry* currentContactGeometry= YADE_CAST<SpheresContactGeometry*>(contact->interactionGeometry.get());
- ViscoelasticInteraction* currentContactPhysics = YADE_CAST<ViscoelasticInteraction*> (contact->interactionPhysics.get());
- if((!currentContactGeometry)||(!currentContactPhysics)) continue;
-
- SimpleViscoelasticBodyParameters* de1 = YADE_CAST<SimpleViscoelasticBodyParameters*>((*bodies)[id1]->physicalParameters.get());
- SimpleViscoelasticBodyParameters* de2 = YADE_CAST<SimpleViscoelasticBodyParameters*>((*bodies)[id2]->physicalParameters.get());
-
- bool isDynamic1 = (*bodies)[id1]->isDynamic;
- bool isDynamic2 = (*bodies)[id2]->isDynamic;
-
- Vector3r& shearForce = currentContactPhysics->shearForce;
-
- if (contact->isNew) shearForce=Vector3r(0,0,0);
-
- Vector3r axis;
- Real angle;
-
- /// Here is the code with approximated rotations ///
-
-
-
- axis = currentContactPhysics->prevNormal.Cross(currentContactGeometry->normal);
- shearForce -= shearForce.Cross(axis);
- //angle = dt*0.5*currentContactGeometry->normal.Dot(de1->angularVelocity+de2->angularVelocity);
- //FIXME: if one body is kinematic then assumed its rotation centre does not lies along the normal
- //(i.e. virtual sphere, which replaces this kinematic body on contact, does not rotate)
- Vector3r summaryAngularVelocity(0,0,0);
- if (isDynamic1) summaryAngularVelocity += de1->angularVelocity;
- if (isDynamic2) summaryAngularVelocity += de2->angularVelocity;
- angle = dt*0.5*currentContactGeometry->normal.Dot(summaryAngularVelocity);
- axis = angle*currentContactGeometry->normal;
- shearForce -= shearForce.Cross(axis);
-
- /// Here is the code with exact rotations ///
-
- // Quaternionr q;
- //
- // axis = currentContactPhysics->prevNormal.cross(currentContactGeometry->normal);
- // angle = acos(currentContactGeometry->normal.dot(currentContactPhysics->prevNormal));
- // q.fromAngleAxis(angle,axis);
- //
- // currentContactPhysics->shearForce = currentContactPhysics->shearForce*q;
- //
- // angle = dt*0.5*currentContactGeometry->normal.dot(de1->angularVelocity+de2->angularVelocity);
- // axis = currentContactGeometry->normal;
- // q.fromAngleAxis(angle,axis);
- // currentContactPhysics->shearForce = q*currentContactPhysics->shearForce;
-
- /// ///
-
- Vector3r x = currentContactGeometry->contactPoint;
- Vector3r c1x = (x - de1->se3.position);
- Vector3r c2x = (x - de2->se3.position);
- /// The following definition of c1x and c2x is to avoid "granular ratcheting"
- /// (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN,
- /// "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena",
- /// ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors)
- Vector3r _c1x_ = (isDynamic1) ? currentContactGeometry->radius1*currentContactGeometry->normal : x - de1->zeroPoint;
- Vector3r _c2x_ = (isDynamic2) ? -currentContactGeometry->radius2*currentContactGeometry->normal : x - de2->zeroPoint;
- Vector3r relativeVelocity = (de2->velocity+de2->angularVelocity.Cross(_c2x_)) - (de1->velocity+de1->angularVelocity.Cross(_c1x_));
- Real normalVelocity = currentContactGeometry->normal.Dot(relativeVelocity);
- Vector3r shearVelocity = relativeVelocity-normalVelocity*currentContactGeometry->normal;
- //Vector3r shearDisplacement = shearVelocity*dt;
- //shearForce -= currentContactPhysics->ks*shearDisplacement;
- shearForce -= (currentContactPhysics->ks*dt+currentContactPhysics->cs)*shearVelocity;
-
- Real un=currentContactGeometry->penetrationDepth;
- //currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal;
- currentContactPhysics->normalForce = ( currentContactPhysics->kn * std::max(un,(Real) 0) - currentContactPhysics->cn * normalVelocity ) * currentContactGeometry->normal;
-
- // PFC3d SlipModel, is using friction angle. CoulombCriterion
- Real maxFs = currentContactPhysics->normalForce.SquaredLength() * std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
- if( shearForce.SquaredLength() > maxFs )
- {
- maxFs = Mathr::Sqrt(maxFs) / shearForce.Length();
- shearForce *= maxFs;
- }
- ////////// PFC3d SlipModel
-
- Vector3r f = currentContactPhysics->normalForce + shearForce;
-
- // it will be some macro( body->physicalActions, ActionType , bodyId )
- static_cast<Force*> ( ncb->physicalActions->find( id1 , actionForceIndex).get() )->force -= f;
- static_cast<Force*> ( ncb->physicalActions->find( id2 , actionForceIndex ).get() )->force += f;
-
- static_cast<Momentum*>( ncb->physicalActions->find( id1 , actionMomentumIndex ).get() )->momentum -= c1x.Cross(f);
- static_cast<Momentum*>( ncb->physicalActions->find( id2 , actionMomentumIndex ).get() )->momentum += c2x.Cross(f);
-
- currentContactPhysics->prevNormal = currentContactGeometry->normal;
- }
- }
-}
-
-YADE_PLUGIN();
Deleted: trunk/pkg/dem/Engine/StandAloneEngine/SimpleViscoelasticContactLaw.hpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/SimpleViscoelasticContactLaw.hpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/dem/Engine/StandAloneEngine/SimpleViscoelasticContactLaw.hpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -1,44 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2004 by Olivier Galizzi *
-* olivier.galizzi@xxxxxxx *
-* *
-* Copyright (C) 2008 by Sergei Dorofeenko *
-* sega@xxxxxxxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/core/InteractionSolver.hpp>
-
-class PhysicalAction;
-
-/// This class provides linear viscoelastic contact model
-class SimpleViscoelasticContactLaw : public InteractionSolver
-{
-/// Attributes
- private :
- shared_ptr<PhysicalAction> actionForce;
- shared_ptr<PhysicalAction> actionMomentum;
- int actionForceIndex;
- int actionMomentumIndex;
-
- public :
-// int sdecGroupMask;
- bool momentRotationLaw;
-
- SimpleViscoelasticContactLaw();
- void action(MetaBody*);
-
- protected :
- void registerAttributes();
- NEEDS_BEX("Force","Momentum");
- REGISTER_CLASS_NAME(SimpleViscoelasticContactLaw);
- REGISTER_BASE_CLASS_NAME(InteractionSolver);
-};
-
-REGISTER_SERIALIZABLE(SimpleViscoelasticContactLaw);
-
-
Modified: trunk/pkg/dem/PreProcessor/MembraneTest.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/MembraneTest.cpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/dem/PreProcessor/MembraneTest.cpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -12,8 +12,9 @@
#include "MembraneTest.hpp"
//#include<yade/pkg-dem/ElasticContactLaw.hpp>
-#include<yade/pkg-dem/SimpleViscoelasticContactLaw.hpp>
+#include<yade/pkg-dem/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp>
//#include<yade/pkg-dem/ElasticBodyParameters2BcpConnection4ElasticContactInteraction.hpp>
+#include<yade/pkg-common/ConstitutiveLawDispatcher.hpp>
#include<yade/pkg-common/ParticleParameters.hpp>
#include<yade/pkg-common/Sphere.hpp>
#include<yade/pkg-common/ef2_BshTube_BssSweptSphereLineSegment_makeBssSweptSphereLineSegment.hpp>
@@ -335,13 +336,16 @@
shared_ptr<PhysicalParametersMetaEngine> orientationIntegrator(new PhysicalParametersMetaEngine);
orientationIntegrator->add("LeapFrogOrientationIntegrator");
+ shared_ptr<ConstitutiveLawDispatcher> constitutiveLaw(new ConstitutiveLawDispatcher);
+ constitutiveLaw->add("ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw");
+
rootBody->engines.clear();
rootBody->engines.push_back(shared_ptr<Engine>(new PhysicalActionContainerReseter));
rootBody->engines.push_back(boundingVolumeDispatcher);
rootBody->engines.push_back(shared_ptr<Engine>(new PersistentSAPCollider));
rootBody->engines.push_back(interactionGeometryDispatcher);
rootBody->engines.push_back(interactionPhysicsDispatcher);
- rootBody->engines.push_back(shared_ptr<Engine>(new SimpleViscoelasticContactLaw));
+ rootBody->engines.push_back(constitutiveLaw);
rootBody->engines.push_back(gravityCondition);
rootBody->engines.push_back(applyActionDispatcher);
rootBody->engines.push_back(positionIntegrator);
Modified: trunk/pkg/dem/PreProcessor/TestSimpleViscoelastic.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/TestSimpleViscoelastic.cpp 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/dem/PreProcessor/TestSimpleViscoelastic.cpp 2009-03-12 17:27:35 UTC (rev 1718)
@@ -30,10 +30,11 @@
#include<yade/pkg-common/PhysicalParametersMetaEngine.hpp>
#include<yade/pkg-common/Sphere.hpp>
#include<yade/pkg-common/Box.hpp>
+#include<yade/pkg-common/ConstitutiveLawDispatcher.hpp>
#include<yade/pkg-dem/RigidBodyRecorder.hpp>
#include<yade/pkg-dem/SimpleViscoelasticSpheresInteractionRecorder.hpp>
#include<yade/pkg-dem/SimpleViscoelasticBodyParameters.hpp>
-#include<yade/pkg-dem/SimpleViscoelasticContactLaw.hpp>
+#include<yade/pkg-dem/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp>
#include<yade/pkg-dem/SimpleViscoelasticRelationships.hpp>
#include<yade/pkg-common/GravityEngines.hpp>
@@ -190,13 +191,16 @@
shared_ptr<PhysicalParametersMetaEngine> orientationIntegrator(new PhysicalParametersMetaEngine);
orientationIntegrator->add("LeapFrogOrientationIntegrator");
+ shared_ptr<ConstitutiveLawDispatcher> constitutiveLaw(new ConstitutiveLawDispatcher);
+ constitutiveLaw->add("ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw");
+
rootBody->engines.clear();
rootBody->engines.push_back(shared_ptr<Engine>(new PhysicalActionContainerReseter));
rootBody->engines.push_back(boundingVolumeDispatcher);
rootBody->engines.push_back(shared_ptr<Engine>(new PersistentSAPCollider));
rootBody->engines.push_back(interactionGeometryDispatcher);
rootBody->engines.push_back(interactionPhysicsDispatcher);
- rootBody->engines.push_back(shared_ptr<Engine>(new SimpleViscoelasticContactLaw));
+ rootBody->engines.push_back(constitutiveLaw);
rootBody->engines.push_back(gravityCondition);
rootBody->engines.push_back(applyActionDispatcher);
rootBody->engines.push_back(positionIntegrator);
Modified: trunk/pkg/dem/SConscript
===================================================================
--- trunk/pkg/dem/SConscript 2009-03-10 13:22:50 UTC (rev 1717)
+++ trunk/pkg/dem/SConscript 2009-03-12 17:27:35 UTC (rev 1718)
@@ -664,7 +664,8 @@
#'ElasticBodyParameters2BcpConnection4ElasticContactInteraction',
'BssSweptSphereLineSegment',
'SimpleViscoelasticRelationships',
- 'SimpleViscoelasticContactLaw',
+ 'ConstitutiveLawDispatcher',
+ 'ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw',
'SimpleViscoelasticBodyParameters',
'ViscoelasticInteraction',
'InteractingSphere','InteractingNode',
@@ -1037,15 +1038,6 @@
,'SpheresContactGeometry'
])
- ,env.SharedLibrary('SimpleViscoelasticContactLaw'
- ,['Engine/StandAloneEngine/SimpleViscoelasticContactLaw.cpp']
- ,LIBS=env['LIBS']+['ViscoelasticInteraction'
- ,'SimpleViscoelasticBodyParameters'
- ,'SpheresContactGeometry'
- ,'Force'
- ,'Momentum'
- ])
-
,env.SharedLibrary('ContactLaw1',
['Engine/StandAloneEngine/ContactLaw1.cpp'],
LIBS=env['LIBS']+['SDECLinkPhysics',
@@ -1078,7 +1070,8 @@
,'Shop'
,'InteractingBox'
,'SimpleViscoelasticRelationships'
- ,'SimpleViscoelasticContactLaw'
+ ,'ConstitutiveLawDispatcher'
+ ,'ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw'
,'SimpleViscoelasticBodyParameters'
,'ViscoelasticInteraction'
,'SimpleViscoelasticSpheresInteractionRecorder'
@@ -1131,8 +1124,8 @@
['Engine/EngineUnit/ef2_Spheres_NormalShear_ElasticFrictionalLaw.cpp'],
LIBS=env['LIBS']+['ConstitutiveLawDispatcher','NormalShearInteractions','SpheresContactGeometry']),
- env.SharedLibrary('Spheres_Viscoelastic_SimpleViscoelasticContactLaw',
- ['Engine/EngineUnit/Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp'],
+ env.SharedLibrary('ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw',
+ ['Engine/EngineUnit/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp'],
LIBS=env['LIBS']+['ConstitutiveLawDispatcher','ViscoelasticInteraction','SpheresContactGeometry','RigidBodyParameters']),
#env.SharedLibrary('MicroMacroAnalyser',
