yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #03224
[Branch ~yade-dev/yade/trunk] Rev 1995: Beginning of adapting code about simple shear (to train with bzr almost)
------------------------------------------------------------
revno: 1995
committer: jduriez <jduriez@c1solimara-l>
branch nick: trunk
timestamp: Mon 2010-02-01 16:20:54 +0100
message:
Beginning of adapting code about simple shear (to train with bzr almost)
- There are still indeed the YADE_REQUIRE_FEATURE(PHYS_PAR) lines which prevent these files to be compiled.
- There is only one change that will be detected by your news compilings : change of a line in ElasticContactLaw. I corrected a link towards an URL corresponding to "old
wiki"
I made "bzr up" just before, so I hope I did nothing bad... :-~
modified:
Yade.kdevelop
pkg/common/Engine/PartialEngine/CinemCNCEngine.cpp
pkg/common/Engine/PartialEngine/CinemCNCEngine.hpp
pkg/common/Engine/PartialEngine/CinemDNCEngine.cpp
pkg/common/Engine/PartialEngine/CinemDNCEngine.hpp
pkg/dem/Engine/Functor/CL1Relationships.cpp
pkg/dem/Engine/GlobalEngine/ContactLaw1.cpp
pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp
pkg/dem/PreProcessor/SimpleShear.cpp
pkg/dem/PreProcessor/SimpleShear.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 'Yade.kdevelop'
--- Yade.kdevelop 2008-11-21 11:05:06 +0000
+++ Yade.kdevelop 2010-02-01 15:20:54 +0000
@@ -10,16 +10,16 @@
<projectname>Yade</projectname>
<projectdirectory>.</projectdirectory>
<absoluteprojectpath>false</absoluteprojectpath>
- <description></description>
- <defaultencoding></defaultencoding>
+ <description/>
+ <defaultencoding/>
<versioncontrol/>
</general>
<kdevcustomproject>
<run>
<mainprogram>/home/vincent/yade-local/trunk</mainprogram>
<directoryradio>executable</directoryradio>
- <programargs></programargs>
- <globaldebugarguments></globaldebugarguments>
+ <programargs/>
+ <globaldebugarguments/>
<globalcwd>/tmp</globalcwd>
<useglobalprogram>false</useglobalprogram>
<terminal>false</terminal>
@@ -30,7 +30,7 @@
</run>
<build>
<buildtool>make</buildtool>
- <builddir></builddir>
+ <builddir/>
</build>
<make>
<abortonerror>false</abortonerror>
@@ -38,8 +38,8 @@
<prio>0</prio>
<dontact>false</dontact>
<makebin>scons</makebin>
- <defaulttarget></defaulttarget>
- <makeoptions></makeoptions>
+ <defaulttarget/>
+ <makeoptions/>
<selectedenvironment>default</selectedenvironment>
<environments>
<default>
@@ -66,9 +66,9 @@
<blacklist/>
<other>
<prio>0</prio>
- <otherbin></otherbin>
- <defaulttarget></defaulttarget>
- <otheroptions></otheroptions>
+ <otherbin/>
+ <defaulttarget/>
+ <otheroptions/>
<selectedenvironment>default</selectedenvironment>
<environments>
<default/>
@@ -77,11 +77,11 @@
</kdevcustomproject>
<kdevdebugger>
<general>
- <dbgshell></dbgshell>
- <gdbpath></gdbpath>
- <configGdbScript></configGdbScript>
- <runShellScript></runShellScript>
- <runGdbScript></runGdbScript>
+ <dbgshell/>
+ <gdbpath/>
+ <configGdbScript/>
+ <runShellScript/>
+ <runGdbScript/>
<breakonloadinglibs>true</breakonloadinglibs>
<separatetty>false</separatetty>
<floatingtoolbar>false</floatingtoolbar>
@@ -154,7 +154,7 @@
<root>/usr/share/qt3</root>
<designerintegration>EmbeddedKDevDesigner</designerintegration>
<qmake>/usr/bin/qmake-qt3</qmake>
- <designer>/usr/bin/designer</designer>
+ <designer></designer>
<designerpluginpaths/>
</qt>
<codecompletion>
@@ -183,7 +183,7 @@
<resolveIncludePathsUsingMakeExperimental>false</resolveIncludePathsUsingMakeExperimental>
</codecompletion>
<creategettersetter>
- <prefixGet></prefixGet>
+ <prefixGet/>
<prefixSet>set</prefixSet>
<prefixVariable>m_,_</prefixVariable>
<parameterName>theValue</parameterName>
=== modified file 'pkg/common/Engine/PartialEngine/CinemCNCEngine.cpp'
--- pkg/common/Engine/PartialEngine/CinemCNCEngine.cpp 2010-01-10 09:09:32 +0000
+++ pkg/common/Engine/PartialEngine/CinemCNCEngine.cpp 2010-02-01 15:20:54 +0000
@@ -8,7 +8,8 @@
#include "CinemCNCEngine.hpp"
-#include<yade/pkg-common/RigidBodyParameters.hpp>
+// #include<yade/pkg-common/RigidBodyParameters.hpp> , remplace par :
+#include<yade/core/State.hpp>
#include<yade/pkg-common/Box.hpp>
#include<yade/pkg-dem/FrictPhys.hpp>
#include<yade/core/Scene.hpp>
@@ -22,6 +23,9 @@
firstRun=true;
shearSpeed=0;
alpha=Mathr::PI/2.0;;
+ gamma_save.resize(0);
+ temoin_save.resize(0);
+ temoin=0;
gamma=0;
gammalim=0;
id_boxhaut=3;
@@ -34,16 +38,15 @@
F_0=0;
Key="";
it_depart=0;
- k=1;
LOG=0;
wallDamping = 0.2;
coeff_dech=1.0;
}
-void CinemCNCEngine::applyCondition(Body * body)
+void CinemCNCEngine::applyCondition(Scene * ncb)
{
- if(LOG) cout << "debut applyCondi !!" << endl;
+ if(LOG) cout << "debut applyCondi du CNCEngine !!" << endl;
leftbox = Body::byId(id_boxleft);
rightbox = Body::byId(id_boxright);
frontbox = Body::byId(id_boxfront);
@@ -51,16 +54,18 @@
topbox = Body::byId(id_boxhaut);
boxbas = Body::byId(id_boxbas);
-
+ if(LOG) cout << "gamma = " << lexical_cast<string>(gamma) << " et gammalim = " << lexical_cast<string>(gammalim) << endl;
if(gamma<=gammalim)
{
- if(temoin==0 || ( (Omega::instance().getCurrentIteration() - it_depart) % k ==0) )
+ if(LOG) cout << "Je suis bien dans la partie gamma < gammalim" << endl;
+ if(temoin==0)
{
if(LOG) cout << "Je veux maintenir la Force a F_0 = : " << F_0 << endl;
- letMove(body);
temoin=1;
}
+ letMove(ncb);
+
}
else if (temoin<2)
{
@@ -75,11 +80,20 @@
Omega::instance().stopSimulationLoop();
}
+ for(unsigned int j=0;j<gamma_save.size();j++)
+ {
+ if ((gamma > gamma_save[j]) && (temoin_save[j]==0))
+ {
+ stopMovement(); // reset of all the speeds before the save
+ Omega::instance().saveSimulation(Key+"_"+lexical_cast<string> (floor(gamma*1000)) +"_" +lexical_cast<string> (floor(gamma*10000)-10*floor(gamma*1000))+ "mmsheared.xml");
+ temoin_save[j]=1;
+ }
+ }
+
}
-void CinemCNCEngine::letMove(Body * body)
+void CinemCNCEngine::letMove(Scene * ncb)
{
- Scene * ncb = YADE_CAST<Scene*>(body);
shared_ptr<BodyContainer> bodies = ncb->bodies;
if(LOG) cout << "It : " << Omega::instance().getCurrentIteration() << endl;
@@ -89,32 +103,27 @@
Real dx = shearSpeed * dt;
- Real Ysup = (topbox->physicalParameters.get())->se3.position.Y();
- Real Ylat = (leftbox->physicalParameters.get())->se3.position.Y();
+ Real Ysup = topbox->state->pos.Y();
+ Real Ylat = leftbox->state->pos.Y();
// Changes in vertical and horizontal position :
- (topbox->physicalParameters.get())->se3.position += Vector3r(dx,deltaU,0);
+ topbox->state->pos += Vector3r(dx,deltaU,0);
- (leftbox->physicalParameters.get())->se3.position += Vector3r(dx/2.0,deltaU/2.0,0);
- (rightbox->physicalParameters.get())->se3.position += Vector3r(dx/2.0,deltaU/2.0,0);
+ leftbox->state->pos += Vector3r(dx/2.0,deltaU/2.0,0);
+ rightbox->state->pos += Vector3r(dx/2.0,deltaU/2.0,0);
if(LOG) cout << "deltaU reellemt applique :" << deltaU << endl;
- if(LOG) cout << "qui nous a emmene en : y = " << ((topbox->physicalParameters.get())->se3.position).Y() << endl;
- if(LOG) cout << "soit un decalage par rapport a position intiale : " << ((topbox->physicalParameters.get())->se3.position.Y()) - Y0 << endl;
+ if(LOG) cout << "qui nous a emmene en : y = " <<(topbox->state->pos).Y() << endl;
+ if(LOG) cout << "soit un decalage par rapport a position intiale : " << (topbox->state->pos.Y()) - Y0 << endl;
- Real Ysup_mod = (topbox->physicalParameters.get())->se3.position.Y();
- Real Ylat_mod = (leftbox->physicalParameters.get())->se3.position.Y();
+ Real Ysup_mod = topbox->state->pos.Y();
+ Real Ylat_mod = leftbox->state->pos.Y();
// with the corresponding velocities :
- RigidBodyParameters * rb = dynamic_cast<RigidBodyParameters*>(topbox->physicalParameters.get());
- rb->velocity = Vector3r(shearSpeed,deltaU/dt,0);
-
- rb = dynamic_cast<RigidBodyParameters*>(leftbox->physicalParameters.get());
- rb->velocity = Vector3r(shearSpeed/2.0,deltaU/(2.0*dt),0);
-
- rb = dynamic_cast<RigidBodyParameters*>(rightbox->physicalParameters.get());
- rb->velocity = Vector3r(shearSpeed/2.0,deltaU/(2.0*dt),0);
+ topbox->state->vel = Vector3r(shearSpeed,deltaU/dt,0);
+ leftbox->state->vel = Vector3r(shearSpeed/2.0,deltaU/(2.0*dt),0);
+ rightbox->state->vel = Vector3r(shearSpeed/2.0,deltaU/(2.0*dt),0);
// Then computation of the angle of the rotation to be done :
computeAlpha();
@@ -132,13 +141,11 @@
qcorr.FromAxisAngle(Vector3r(0,0,1),dalpha);
// On applique la rotation en changeant l'orientation des plaques, leurs vang et en affectant donc alpha
- rb = dynamic_cast<RigidBodyParameters*>(leftbox->physicalParameters.get());
- rb->se3.orientation = qcorr*rb->se3.orientation;
- rb->angularVelocity = Vector3r(0,0,1)*dalpha/dt;
+ leftbox->state->ori = qcorr*leftbox->state->ori;
+ leftbox->state->angVel = Vector3r(0,0,1)*dalpha/dt;
- rb = dynamic_cast<RigidBodyParameters*>(rightbox->physicalParameters.get());
- rb->se3.orientation = qcorr*rb->se3.orientation;
- rb->angularVelocity = Vector3r(0,0,1)*dalpha/dt;
+ rightbox->state->ori = qcorr*rightbox->state->ori;
+ rightbox->state->angVel = Vector3r(0,0,1)*dalpha/dt;
gamma+=dx;
}
@@ -147,8 +154,8 @@
void CinemCNCEngine::computeAlpha()
{
Quaternionr orientationLeftBox,orientationRightBox;
- orientationLeftBox = (dynamic_cast<RigidBodyParameters*>(leftbox->physicalParameters.get()) )->se3.orientation;
- orientationRightBox = (dynamic_cast<RigidBodyParameters*>(rightbox->physicalParameters.get()) )->se3.orientation;
+ orientationLeftBox = leftbox->state->ori;
+ orientationRightBox = rightbox->state->ori;
if(orientationLeftBox!=orientationRightBox)
{
cout << "WARNING !!! your lateral boxes have not the same orientation, you're not in the case of a box imagined for creating these engines" << endl;
@@ -183,7 +190,7 @@
it_depart = Omega::instance().getCurrentIteration();
alpha=Mathr::PI/2.0;;
- Y0 = (topbox->physicalParameters.get())->se3.position.Y();
+ Y0 = topbox->state->pos.Y();
cout << "Y0 initialise à : " << Y0 << endl;
F_0 = F_sup.Y();
cout << "F_0 initialise à : " << F_0 << endl;
@@ -200,7 +207,7 @@
}
else
{
- Real Ycourant = (topbox->physicalParameters.get())->se3.position.Y();
+ Real Ycourant = topbox->state->pos.Y();
deltaU = ( F_sup.Y() - F_0 )/(stiffness);
if(LOG) cout << "Lors du calcul de DU (utile pour deltaU) : F_0 = " << F_0 << "; Y0 = " << Y0 << "; Ycourant = " << Ycourant << endl;
}
@@ -235,18 +242,15 @@
void CinemCNCEngine::stopMovement()
{
// annulation de la vitesse de la plaque du haut
- RigidBodyParameters * rb = YADE_CAST<RigidBodyParameters*>(topbox->physicalParameters.get());
- rb->velocity = Vector3r(0,0,0);
+ topbox->state->vel = Vector3r(0,0,0);
// de la plaque gauche
- rb = YADE_CAST<RigidBodyParameters*>(leftbox->physicalParameters.get());
- rb->velocity = Vector3r(0,0,0);
- rb->angularVelocity = Vector3r(0,0,0);
+ leftbox->state->vel = Vector3r(0,0,0);
+ leftbox->state->angVel = Vector3r(0,0,0);
// de la plaque droite
- rb = YADE_CAST<RigidBodyParameters*>(rightbox->physicalParameters.get());
- rb->velocity = Vector3r(0,0,0);
- rb->angularVelocity = Vector3r(0,0,0);
+ rightbox->state->vel = Vector3r(0,0,0);
+ rightbox->state->angVel = Vector3r(0,0,0);
}
void CinemCNCEngine::computeStiffness(Scene* ncb)
@@ -283,5 +287,4 @@
-YADE_REQUIRE_FEATURE(PHYSPAR);
=== modified file 'pkg/common/Engine/PartialEngine/CinemCNCEngine.hpp'
--- pkg/common/Engine/PartialEngine/CinemCNCEngine.hpp 2009-12-09 17:11:51 +0000
+++ pkg/common/Engine/PartialEngine/CinemCNCEngine.hpp 2010-02-01 15:20:54 +0000
@@ -47,7 +47,6 @@
int temoin,it_stop
,it_depart // the number of the first it at which plates are moved
- ,k // plates are moved only all "k" iterations to let the sample relax
;
Vector3r prevF_sup; // the force acting on the upper plate at the previous time step
@@ -61,7 +60,7 @@
shared_ptr<Body> boxbas;
public :
CinemCNCEngine();
- void applyCondition(Body * body)
+ void applyCondition(Scene * body)
,computeAlpha()
;
@@ -71,6 +70,10 @@
,wallDamping
;
+ std::vector<Real> gamma_save // the values of gamma, in meters, at which a save is performed
+ ,temoin_save
+ ;
+
bool LOG; //controls messages output on screen
body_id_t id_boxhaut, // the id of the upper wall : defined in the constructor
@@ -84,8 +87,8 @@
protected :
- REGISTER_ATTRIBUTES(PartialEngine,(shearSpeed)(gammalim)(prevF_sup)(firstRun)(id_boxhaut)(id_boxbas)(id_boxleft)(id_boxright)(id_boxfront)(id_boxback)(Y0)(F_0)(k)(max_vel)(wallDamping)(Key)(LOG)(coeff_dech));
- void letMove(Body* body);
+ REGISTER_ATTRIBUTES(PartialEngine,(shearSpeed)(gammalim)(gamma)(gamma_save)(temoin_save)(prevF_sup)(firstRun)(id_boxhaut)(id_boxbas)(id_boxleft)(id_boxright)(id_boxfront)(id_boxback)(Y0)(F_0)(max_vel)(wallDamping)(Key)(LOG)(coeff_dech));
+ void letMove(Scene* body);
void computeDu(Scene* ncb);
void stopMovement(); // to cancel all the velocities when gammalim is reached
void computeStiffness(Scene* ncb);
=== modified file 'pkg/common/Engine/PartialEngine/CinemDNCEngine.cpp'
--- pkg/common/Engine/PartialEngine/CinemDNCEngine.cpp 2009-12-04 23:07:34 +0000
+++ pkg/common/Engine/PartialEngine/CinemDNCEngine.cpp 2010-02-01 15:20:54 +0000
@@ -14,89 +14,139 @@
#include"CinemDNCEngine.hpp"
-CinemDNCEngine::CinemDNCEngine()
+CinemDNCEngine::CinemDNCEngine() : leftbox(new Body), rightbox(new Body), topbox(new Body)
{
- rotationAxis=Vector3r(0,0,1);
+ gamma_save.resize(0);
+ temoin_save.resize(0);
+ temoinfin=0;
shearSpeed=0;
- theta=0;
- thetalim=9.0/10.0*Mathr::PI/2.0;
+ gamma=0;
+ gammalim=0;
id_boxhaut=3;
+ id_boxleft=0;
+ id_boxright=2;
}
void CinemDNCEngine::applyCondition(Scene * body)
{
- if(theta<=thetalim)
- {
- applyRotTranslation(body); // to let move the lateral walls
- applyTranslation(body); // to let move the upper wall
+ leftbox = Body::byId(id_boxleft);
+ rightbox = Body::byId(id_boxright);
+ topbox = Body::byId(id_boxhaut);
+
+ if( ((shearSpeed > 0) && (gamma<=gammalim)) || ((shearSpeed < 0) /*&& (gamma>=gammalim)*/ ) )
+ {
+ if(temoinfin!=0)
+ temoinfin=0;
+ letMove(ncb);
+ }
+ else
+ {
+ stopMovement();
+ if(temoinfin==0)
+ {
+ Omega::instance().saveSimulation(Key + "finCis.xml");
+ temoinfin=1;
+ }
+ }
+
+ for(unsigned int j=0;j<gamma_save.size();j++)
+ {
+ if ( ( ( (shearSpeed>0)&&(gamma > gamma_save[j]) ) || ((shearSpeed<0)&&(gamma < gamma_save[j])) ) && (temoin_save[j]==0) )
+ {
+ stopMovement(); // reset of all the speeds before the save
+ Omega::instance().saveSimulation(Key+"_"+lexical_cast<string> (floor(gamma*1000)) + "mmsheared.xml");
+ temoin_save[j]=1;
+ }
}
}
-
-void CinemDNCEngine::applyRotTranslation(Scene * ncb)
+void CinemDNCEngine::letMove(MetaBody * ncb)
{
shared_ptr<BodyContainer> bodies = ncb->bodies;
-
- Yplaqsup=((*bodies)[id_boxhaut]->physicalParameters.get())->se3.position.Y(); // the height of the sample, which may be different from the initial "height" defined in the Preprocessor
+ Real dt = Omega::instance().getTimeStep();
+ Real dx = shearSpeed * dt;
+
+ (topbox->physicalParameters.get())->se3.position += Vector3r(dx,0,0);
+
+ (leftbox->physicalParameters.get())->se3.position += Vector3r(dx/2.0,0,0);
+ (rightbox->physicalParameters.get())->se3.position += Vector3r(dx/2.0,0,0);
+
+ Real Ysup = (topbox->physicalParameters.get())->se3.position.Y();
+ Real Ylat = (leftbox->physicalParameters.get())->se3.position.Y();
+
+
+// with the corresponding velocities :
+ RigidBodyParameters * rb = dynamic_cast<RigidBodyParameters*>(topbox->physicalParameters.get());
+ rb->velocity = Vector3r(shearSpeed,0,0);
+
+ rb = dynamic_cast<RigidBodyParameters*>(leftbox->physicalParameters.get());
+ rb->velocity = Vector3r(shearSpeed/2.0,0,0);
+
+ rb = dynamic_cast<RigidBodyParameters*>(rightbox->physicalParameters.get());
+ rb->velocity = Vector3r(shearSpeed/2.0,0,0);
+
+// Then computation of the angle of the rotation to be done :
+ computeAlpha();
+ if (alpha == Mathr::PI/2.0) // Case of the very beginning
+ {
+ dalpha = - Mathr::ATan( dx / (Ysup -Ylat) );
+ }
+ else
+ {
+ Real A = (Ysup - Ylat) * 2.0*Mathr::Tan(alpha) / (2.0*(Ysup - Ylat) + dx*Mathr::Tan(alpha) );
+ dalpha = Mathr::ATan( (A - Mathr::Tan(alpha))/(1.0 + A * Mathr::Tan(alpha)));
+ }
- subscribedBodies.clear();
- subscribedBodies.push_back(0); // 0 and 2 are in my case the ids of the lateral walls
- subscribedBodies.push_back(2);
- std::vector<int>::const_iterator ii = subscribedBodies.begin();
- std::vector<int>::const_iterator iiEnd = subscribedBodies.end();
-
- Real dt = Omega::instance().getTimeStep();
- Real dx = shearSpeed*dt;
- Real dtheta;
- dtheta=Mathr::ATan(dx/(Yplaqsup+dx*Mathr::Tan(theta)+Yplaqsup*Mathr::Pow(Mathr::Tan(theta),2)));
- Quaternionr q;
- q.FromAxisAngle(rotationAxis,-dtheta);
-
- for(;ii!=iiEnd;++ii)
- {
- RigidBodyParameters * rb = dynamic_cast<RigidBodyParameters*>((*bodies)[*ii]->physicalParameters.get());
-
-// La partie de rotation :
- rb->se3.orientation = q*rb->se3.orientation;
- rb->se3.orientation.Normalize();
-/* rb->angularVelocity = rotationAxis*angularVelocity;*/
- rb->velocity = Vector3r(0,0,0);
-// La partie de translation :
- Real YcentreW1; //the altitude of the center of both lateral walls
- YcentreW1=((*bodies)[0]->physicalParameters.get())->se3.position.Y();
- rb->se3.position+=dt*(shearSpeed*YcentreW1/Yplaqsup)*Vector3r(1,0,0);
- }
- theta+=dtheta;
-}
-
-
-void CinemDNCEngine::applyTranslation(Scene * ncb)
-{
- shared_ptr<BodyContainer>& bodies = ncb->bodies;
-
- int id_boxhaut=3;
- Real dt = Omega::instance().getTimeStep();
-
- if(ParticleParameters* p = dynamic_cast<ParticleParameters*>((*bodies)[id_boxhaut]->physicalParameters.get()))
- {
- p->se3.position += dt*shearSpeed*Vector3r(1,0,0); // by def of shearSpeed
- p->velocity = shearSpeed*Vector3r(1,0,0);
- }
- else if(PhysicalParameters* b = dynamic_cast<PhysicalParameters*>((*bodies)[id_boxhaut]->physicalParameters.get()))
- { // NOT everyone has velocity !
- b->se3.position += dt*shearSpeed*Vector3r(1,0,0);
- }
- else
- {
- std::cerr << "TranslationEngine::applyCondition, WARNING! dynamic_cast failed! for id: " << id_boxhaut << std::endl;
- }
-
-}
-
-
+ Quaternionr qcorr;
+ qcorr.FromAxisAngle(Vector3r(0,0,1),dalpha);
+
+// On applique la rotation en changeant l'orientation des plaques, leurs vang et en affectant donc alpha
+ rb = dynamic_cast<RigidBodyParameters*>(leftbox->physicalParameters.get());
+ rb->se3.orientation = qcorr*rb->se3.orientation;
+ rb->angularVelocity = Vector3r(0,0,1)*dalpha/dt;
+
+ rb = dynamic_cast<RigidBodyParameters*>(rightbox->physicalParameters.get());
+ rb->se3.orientation = qcorr*rb->se3.orientation;
+ rb->angularVelocity = Vector3r(0,0,1)*dalpha/dt;
+
+ gamma+=dx;
+}
+
+
+void CinemDNCEngine::computeAlpha()
+{
+ Quaternionr orientationLeftBox,orientationRightBox;
+ orientationLeftBox = (dynamic_cast<RigidBodyParameters*>(leftbox->physicalParameters.get()) )->se3.orientation;
+ orientationRightBox = (dynamic_cast<RigidBodyParameters*>(rightbox->physicalParameters.get()) )->se3.orientation;
+ if(orientationLeftBox!=orientationRightBox)
+ {
+ cout << "WARNING !!! your lateral boxes have not the same orientation, you're not in the case of a box imagined for creating these engines" << endl;
+ }
+ Vector3r axis;
+ Real angle;
+ orientationLeftBox.ToAxisAngle(axis,angle);
+ alpha=Mathr::PI/2.0-angle; // right if the initial orientation of the body (on the beginning of the simulation) is q =(1,0,0,0) = FromAxisAngle((0,0,1),0)
+}
+
+void CinemDNCEngine::stopMovement()
+{
+ // annulation de la vitesse de la plaque du haut
+ RigidBodyParameters * rb = YADE_CAST<RigidBodyParameters*>(topbox->physicalParameters.get());
+ rb->velocity = Vector3r(0,0,0);
+
+ // de la plaque gauche
+ rb = YADE_CAST<RigidBodyParameters*>(leftbox->physicalParameters.get());
+ rb->velocity = Vector3r(0,0,0);
+ rb->angularVelocity = Vector3r(0,0,0);
+
+ // de la plaque droite
+ rb = YADE_CAST<RigidBodyParameters*>(rightbox->physicalParameters.get());
+ rb->velocity = Vector3r(0,0,0);
+ rb->angularVelocity = Vector3r(0,0,0);
+}
YADE_REQUIRE_FEATURE(PHYSPAR);
=== modified file 'pkg/common/Engine/PartialEngine/CinemDNCEngine.hpp'
--- pkg/common/Engine/PartialEngine/CinemDNCEngine.hpp 2009-12-09 17:11:51 +0000
+++ pkg/common/Engine/PartialEngine/CinemDNCEngine.hpp 2010-02-01 15:20:54 +0000
@@ -15,31 +15,45 @@
/*! \brief To apply a zero normal displacement shear for a parallelogram box
-This engine, used in simulations issued from "DirectShearCis" Preprocessor, allows to translate horizontally the upper plate while the lateral ones rotate so that they always keep contact with the lower and upper walls
+This engine, used in simulations issued from "SimpleShear" Preprocessor, allows to translate horizontally the upper plate while the lateral ones rotate so that they always keep contact with the lower and upper walls
*/
class CinemDNCEngine : public PartialEngine
{
private :
- Real Yplaqsup// height of the upper wall, the engine cares to find its value
- ,theta; // the angle between a lateral plate and its original orienation : will increase from 0 to thetalim
+ Real alpha // angle from the lower plate to the left box (trigo wise), the Engine finds itself its value
+ ,dalpha // the increment over alpha, due to vertical displacement of upper box
+ ,gamma // horizontal displacement done since the launch of the calcul first step
+ ;
+
+ int temoinfin;
+ shared_ptr<Body> leftbox;
+ shared_ptr<Body> rightbox;
+ shared_ptr<Body> topbox;
+
+ void letMove(Scene* body);
+ void stopMovement(); // to cancel all the velocities when gammalim is reached
+
public :
CinemDNCEngine();
- void applyCondition(Scene * body);
+ void applyCondition(Scene * body),
+ computeAlpha();
Real shearSpeed // to be defined in the PreProcessor
- ,thetalim // the maximum value of theta, at wich the displacement is stopped
+ ,gammalim // the maximum value of gamma (tangential displacemt), in meters, at wich the displacement is stopped
;
- body_id_t id_boxhaut; // the id of the upper wall : defined in the constructor
- Vector3r rotationAxis; // defined in the constructor
-
+ body_id_t id_boxhaut, // the id of the upper wall : defined in the constructor
+ id_boxleft,
+ id_boxright;
+ std::vector<Real> gamma_save // the values of gamma, in meters, at which a save is performed
+ ,temoin_save
+ ;
+ string Key;
protected :
- void applyRotTranslation(Scene *); // to let move (rotation combined with translation) the lateral walls
- void applyTranslation(Scene *); // to let move (translation) the upper wall
- REGISTER_ATTRIBUTES(PartialEngine,(shearSpeed)(rotationAxis)(theta)(thetalim)(id_boxhaut));
+ REGISTER_ATTRIBUTES(PartialEngine,(shearSpeed)(gammalim)(gamma)(gamma_save)(temoin_save)(id_boxhaut)(id_boxleft)(id_boxright)(Key));
REGISTER_CLASS_NAME(CinemDNCEngine);
REGISTER_BASE_CLASS_NAME(PartialEngine);
};
=== modified file 'pkg/dem/Engine/Functor/CL1Relationships.cpp'
--- pkg/dem/Engine/Functor/CL1Relationships.cpp 2009-12-13 20:11:31 +0000
+++ pkg/dem/Engine/Functor/CL1Relationships.cpp 2010-02-01 15:20:54 +0000
@@ -191,4 +191,3 @@
YADE_PLUGIN((CL1Relationships));
YADE_REQUIRE_FEATURE(PHYSPAR);
-
=== modified file 'pkg/dem/Engine/GlobalEngine/ContactLaw1.cpp'
--- pkg/dem/Engine/GlobalEngine/ContactLaw1.cpp 2009-12-25 14:46:48 +0000
+++ pkg/dem/Engine/GlobalEngine/ContactLaw1.cpp 2010-02-01 15:20:54 +0000
@@ -74,7 +74,7 @@
Real Fn; // la valeur de Fn qui va etre calculee selon différentes manieres puis affectee
Real un = currentContactGeometry->penetrationDepth; // compte positivement lorsq vraie interpenetration
-// cout << "un = " << un << " previousun = " << previousun << " previousFn =" << previousFn << endl;
+// cout << "un = " << un << " alors que unMax = "<< currentContactPhysics->unMax << " et previousun = " << previousun << " et previousFn =" << previousFn << endl;
if(un > currentContactPhysics->unMax) // on est en charge, et au delà et sur la "droite principale"
{
Fn = kn*un;
@@ -219,16 +219,6 @@
///// Real angle_twist(angle * axis.Dot(currentContactGeometry->normal) );
///// Vector3r axis_twist(angle_twist * currentContactGeometry->normal);
//////* no creep
-////////////////////////////// CREEP START ///////////////////////////
-///// Real viscosity_twist = viscosity * std::pow((2 * std::min(currentContactGeometry->radius1,currentContactGeometry->radius2)),2) / 16.0;
-///// Real angle_twist_creeped = angle_twist * (1 - dt/viscosity_twist);
-///// Quaternionr q_twist(currentContactGeometry->normal , angle_twist);
-///// //Quaternionr q_twist_creeped(currentContactGeometry->normal , angle_twist*0.996);
-///// Quaternionr q_twist_creeped(currentContactGeometry->normal , angle_twist_creeped);
-///// Quaternionr q_twist_delta(q_twist_creeped * q_twist.Conjugate() );
-///// // modify the initialRelativeOrientation to substract some twisting
-///// currentContactPhysics->initialRelativeOrientation = currentContactPhysics->initialRelativeOrientation * q_twist_delta;
-////////////////////////////// CREEP END /////////////////////////////
/////*/
///// Vector3r moment_twist(axis_twist * currentContactPhysics->kr);
/////
@@ -304,6 +294,7 @@
}
}
}
+// cout << " En tout :"<< nbreInteracTot << "interactions (reelles)" << endl;
//cerr << "ncount= " << ncount << endl;//REMOVE
// cout << momentAlwaysElastic << endl;
// cout << "Sur " << nbreInteracTot << " interactions (reelles) " << nbreInteracMomPlastif << " se sont vues corriger leur moment" << endl;
=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp 2010-01-21 07:53:17 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp 2010-02-01 15:20:54 +0000
@@ -49,7 +49,7 @@
This class serves also as tutorial and is documented in detail at
- http://yade.wikia.com/wiki/ConstitutiveLawHowto
+ https://yade-dem.org/index.php/ConstitutiveLawHowto
*/
class Law2_Dem3DofGeom_FrictPhys_Basic: public LawFunctor{
public:
=== modified file 'pkg/dem/PreProcessor/SimpleShear.cpp'
--- pkg/dem/PreProcessor/SimpleShear.cpp 2009-12-25 14:46:48 +0000
+++ pkg/dem/PreProcessor/SimpleShear.cpp 2010-02-01 15:20:54 +0000
@@ -21,14 +21,14 @@
#include<yade/pkg-dem/GlobalStiffnessTimeStepper.hpp>
#include <yade/pkg-dem/PositionOrientationRecorder.hpp>
-#include <yade/pkg-dem/PositionRecorder.hpp>
-#include <yade/pkg-dem/PositionSnapshot.hpp>
-#include <yade/pkg-dem/ForceRecorder.hpp>
-#include <yade/pkg-dem/ForceSnapshot.hpp>
+// #include <yade/pkg-dem/PositionRecorder.hpp>
+// #include <yade/pkg-dem/PositionSnapshot.hpp>
+// #include <yade/pkg-dem/ForceRecorder.hpp>
+// #include <yade/pkg-dem/ForceSnapshot.hpp>
-#include<yade/pkg-common/BoxModel.hpp>
+// #include<yade/pkg-common/BoxModel.hpp>
#include<yade/pkg-common/Aabb.hpp>
-#include<yade/pkg-common/SphereModel.hpp>
+// #include<yade/pkg-common/SphereModel.hpp>
#include<yade/core/Scene.hpp>
#include<yade/pkg-common/InsertionSortCollider.hpp>
#include<yade/lib-serialization/IOFormatManager.hpp>
@@ -53,7 +53,7 @@
#include <utility>
using namespace std;
-YADE_REQUIRE_FEATURE(geometricalmodel);
+// YADE_REQUIRE_FEATURE(geometricalmodel);
SimpleShear::SimpleShear () : FileGenerator()
{
@@ -105,10 +105,9 @@
{
rootBody = shared_ptr<Scene>(new Scene);
createActors(rootBody);
- positionRootBody(rootBody);
-
-
-// BoxModel walls
+
+
+// Box walls
shared_ptr<Body> w1; // The left one :
createBox(w1,Vector3r(-thickness/2.0,(height)/2.0,0),Vector3r(thickness/2.0,5*(height/2.0+thickness),profondeur/2.0));
rootBody->bodies->insert(w1);
@@ -116,7 +115,7 @@
shared_ptr<Body> w2; // The lower one :
createBox(w2,Vector3r(width/2.0,-thickness/2.0,0),Vector3r(width/2.0,thickness/2.0,profondeur/2.0));
- YADE_CAST<CohesiveFrictionalMat*>(w2->physicalParameters.get())->frictionAngle = sphereFrictionDeg * Mathr::PI/180.0;; // so that we have phi(spheres-inferior wall)=phi(sphere-sphere)
+ YADE_PTR_CAST<FrictMat> (w2->material)->frictionAngle = sphereFrictionDeg * Mathr::PI/180.0; // so that we have phi(spheres-inferior wall)=phi(sphere-sphere)
rootBody->bodies->insert(w2);
shared_ptr<Body> w3; // The right one
@@ -125,15 +124,14 @@
shared_ptr<Body> w4; // The upper one
createBox(w4,Vector3r(width/2.0,height+thickness/2.0,0),Vector3r(width/2.0,thickness/2.0,profondeur/2.0));
- YADE_CAST<CohesiveFrictionalMat*>(w4->physicalParameters.get())->frictionAngle = sphereFrictionDeg * Mathr::PI/180.0;; // so that we have phi(spheres-superior wall)=phi(sphere-sphere)
- rootBody->bodies->insert(w4);
+ YADE_PTR_CAST<FrictMat> (w4->material)->frictionAngle = sphereFrictionDeg * Mathr::PI/180.0; // so that we have phi(spheres-superior wall)=phi(sphere-sphere)
-// To close the front and the bottom of the box
- shared_ptr<Body> w5;
+// To close the front and the back of the box
+ shared_ptr<Body> w5; // behind
createBox(w5,Vector3r(width/2.0,height/2.0,-profondeur/2.0-thickness/2.0), Vector3r(2.5*width/2.0,height/2.0+thickness,thickness/2.0));
rootBody->bodies->insert(w5);
- shared_ptr<Body> w6;
+ shared_ptr<Body> w6; // the front
createBox(w6,Vector3r(width/2.0,height/2.0,profondeur/2.0+thickness/2.0),
Vector3r(2.5*width/2.0,height/2.0+thickness,thickness/2.0));
rootBody->bodies->insert(w6);
@@ -165,52 +163,56 @@
void SimpleShear::createSphere(shared_ptr<Body>& body, Vector3r position, Real radius)
{
body = shared_ptr<Body>(new Body(0,1));
- shared_ptr<CohesiveFrictionalMat> physics(new CohesiveFrictionalMat);
+ shared_ptr<CohesiveFrictionalMat> mat(new CohesiveFrictionalMat);
shared_ptr<Aabb> aabb(new Aabb);
- shared_ptr<SphereModel> gSphere(new SphereModel);
+// shared_ptr<SphereModel> gSphere(new SphereModel);
shared_ptr<Sphere> iSphere(new Sphere);
- Quaternionr q;
+ Quaternionr q; // to define the initial orientation of the sphere
q.FromAxisAngle( Vector3r(0,0,1),0);
body->isDynamic = true;
-
- physics->angularVelocity = Vector3r(0,0,0);
- physics->velocity = Vector3r(0,0,0);
- physics->mass = 4.0/3.0*Mathr::PI*radius*radius*radius*density;
- physics->inertia = Vector3r(2.0/5.0*physics->mass*radius*radius,2.0/5.0*physics->mass*radius*radius,2.0/5.0*physics->mass*radius*radius);
- physics->se3 = Se3r(position,q);
- physics->young = sphereYoungModulus;
- physics->poisson = spherePoissonRatio;
- physics->frictionAngle = sphereFrictionDeg * Mathr::PI/180.0;
- physics->isCohesive = 1;
+ body->state->pos =position;
+ body->state->ori =q;
+ body->state->vel =Vector3r(0,0,0);
+ body->state->angVel =Vector3r(0,0,0);
+
+ Real masse =4.0/3.0*Mathr::PI*radius*radius*radius*density;
+ body->state->mass =masse;
+ body->state->inertia = Vector3r(2.0/5.0*masse*radius*radius,2.0/5.0*masse*radius*radius,2.0/5.0*masse*radius*radius);
+
+ mat->young = sphereYoungModulus;
+ mat->poisson = spherePoissonRatio;
+ mat->frictionAngle = sphereFrictionDeg * Mathr::PI/180.0;
+ mat->isCohesive = 1;
+ body->material = mat;
aabb->diffuseColor = Vector3r(0,1,0);
- gSphere->radius = radius;
+/* gSphere->radius = radius;
// de quoi avoir des bandes (huit en largeur) de couleur differentes :
gSphere->diffuseColor = ((int)(Mathr::Floor(8*position.X()/width)))%2?Vector3r(0.7,0.7,0.7):Vector3r(0.45,0.45,0.45);
gSphere->wire = false;
- gSphere->shadowCaster = true;
+ gSphere->shadowCaster = true;*/
iSphere->radius = radius;
iSphere->diffuseColor = Vector3r(0.8,0.3,0.3);
- body->shape = iSphere;
- body->geometricalModel = gSphere;
- body->bound = aabb;
- body->physicalParameters = physics;
+ body->shape = iSphere;
+// body->geometricalModel = gSphere;
+ body->bound = aabb;
}
void SimpleShear::createBox(shared_ptr<Body>& body, Vector3r position, Vector3r extents)
{
body = shared_ptr<Body>(new Body(0,1));
- shared_ptr<CohesiveFrictionalMat> physics(new CohesiveFrictionalMat);
+ shared_ptr<CohesiveFrictionalMat> mat(new CohesiveFrictionalMat);
shared_ptr<Aabb> aabb(new Aabb);
- shared_ptr<BoxModel> gBox(new BoxModel);
+// shared_ptr<BoxModel> gBox(new BoxModel);
shared_ptr<Box> iBox(new Box);
+
Quaternionr q;
@@ -218,53 +220,48 @@
body->isDynamic = false;
- physics->angularVelocity = Vector3r(0,0,0);
- physics->velocity = Vector3r(0,0,0);
- physics->mass = extents[0]*extents[1]*extents[2]*density*2;
- physics->inertia = Vector3r(
- physics->mass*(extents[1]*extents[1]+extents[2]*extents[2])/3
- , physics->mass*(extents[0]*extents[0]+extents[2]*extents[2])/3
- , physics->mass*(extents[1]*extents[1]+extents[0]*extents[0])/3
- );
- //physics->mass = 0;
- //physics->inertia = Vector3r(0,0,0);
- physics->se3 = Se3r(position,q);
- physics->young = boxYoungModulus;
- physics->poisson = boxPoissonRatio;
- physics->frictionAngle = 0.0; //default value, modified after for w2 and w4 to have good values of phi(sphere-walls)
- physics->isCohesive = 1;
+ body->state->angVel = Vector3r(0,0,0);
+ body->state->vel = Vector3r(0,0,0);
+// NB : mass and inertia not defined because not used, since Box are not dynamics
+ body->state->pos = position;
+ body->state->ori = q;
+
+ mat->young = boxYoungModulus;
+ mat->poisson = boxPoissonRatio;
+ mat->frictionAngle = 0.0; //default value, modified after for w2 and w4 to have good values of phi(sphere-walls)
+ mat->isCohesive = 1;
+ body->material = mat;
aabb->diffuseColor = Vector3r(1,0,0);
- gBox->extents = extents;
+/* gBox->extents = extents;
gBox->diffuseColor = Vector3r(1,0,0);
gBox->wire = true;
- gBox->shadowCaster = false;
+ gBox->shadowCaster = false;*/
iBox->extents = extents;
iBox->diffuseColor = Vector3r(1,0,0);
- body->bound = aabb;
- body->shape = iBox;
- body->geometricalModel = gBox;
- body->physicalParameters = physics;
+ body->bound = aabb;
+ body->shape = iBox;
+// body->geometricalModel = gBox;
}
void SimpleShear::createActors(shared_ptr<Scene>& rootBody)
{
- shared_ptr<PositionSnapshot> possnap = shared_ptr<PositionSnapshot>(new PositionSnapshot);
- possnap->list_id.clear();
- possnap->list_id.push_back(5);
- possnap->list_id.push_back(10);
- possnap->outputFile="../data/PosSnapshot";
-
- shared_ptr<ForceSnapshot> forcesnap = shared_ptr<ForceSnapshot>(new ForceSnapshot);
- forcesnap->list_id.clear();
- forcesnap->list_id.push_back(5);
- forcesnap->list_id.push_back(10);
- forcesnap->outputFile="../data/ForceSnapshot";
+// shared_ptr<PositionSnapshot> possnap = shared_ptr<PositionSnapshot>(new PositionSnapshot);
+// possnap->list_id.clear();
+// possnap->list_id.push_back(5);
+// possnap->list_id.push_back(10);
+// possnap->outputFile="../data/PosSnapshot";
+//
+// shared_ptr<ForceSnapshot> forcesnap = shared_ptr<ForceSnapshot>(new ForceSnapshot);
+// forcesnap->list_id.clear();
+// forcesnap->list_id.push_back(5);
+// forcesnap->list_id.push_back(10);
+// forcesnap->outputFile="../data/ForceSnapshot";
shared_ptr<CinemDNCEngine> kinematic = shared_ptr<CinemDNCEngine>(new CinemDNCEngine);
kinematic->shearSpeed = shearSpeed;
@@ -314,28 +311,6 @@
}
-void SimpleShear::positionRootBody(shared_ptr<Scene>& rootBody)
-{
- rootBody->isDynamic = false;
-
- Quaternionr q;
- q.FromAxisAngle( Vector3r(0,0,1),0);
-
- shared_ptr<ParticleParameters> physics(new ParticleParameters); // FIXME : fix indexable class PhysicalParameters
- physics->se3 = Se3r(Vector3r(0,0,0),q);
- physics->mass = 0;
- physics->velocity = Vector3r(0,0,0);
- physics->acceleration = Vector3r::ZERO;
-
-
- shared_ptr<Aabb> aabb(new Aabb);
- aabb->diffuseColor = Vector3r(0,0,1);
-
- rootBody->bound = YADE_PTR_CAST<Bound>(aabb);
- rootBody->physicalParameters = physics;
-}
-
-
string SimpleShear::GenerateCloud(vector<BasicSphere>& sphere_list,Vector3r lowerCorner,Vector3r upperCorner,long number,Real rad_std_dev, Real porosity)
{
sphere_list.clear();
@@ -385,18 +360,31 @@
if(importFilename.size() != 0 && filesystem::exists(importFilename) )
{
ifstream loadFile(importFilename.c_str()); // cree l'objet loadFile de la classe ifstream qui va permettre de lire ce qu'il y a dans importFilename
- Real zJF;
+
+// Real zJF;
+// while( !loadFile.eof() ) // tant qu'on n'est pas a la fin du fichier
+// {
+// BasicSphere s; // l'elt de la liste sphere_list (= la sphere) que l'on va lire maintenant
+// loadFile >> s.first.X();// le X de la position de son centre
+// loadFile >> zJF;
+// s.first.Z()=zJF - profondeur/2.0;// le Z de la position de son centre
+// loadFile >> s.first.Y();// le Y de la position de son centre
+// loadFile >> s.second; // son rayon
+// sphere_list.push_back(s);
+// nombre++;
+// }
+ Real it;
while( !loadFile.eof() ) // tant qu'on n'est pas a la fin du fichier
{
BasicSphere s; // l'elt de la liste sphere_list (= la sphere) que l'on va lire maintenant
- loadFile >> s.first.X();// le X de la position de son centre
- loadFile >> zJF;
- s.first.Z()=zJF - profondeur/2.0;// le Z de la position de son centre
+ loadFile >> it;
+ loadFile >> s.second; // son rayon
+ loadFile >> s.first.X();
loadFile >> s.first.Y();// le Y de la position de son centre
- loadFile >> s.second; // son rayon
+ loadFile >> s.first.Z();// le Z de la position de son centre
sphere_list.push_back(s);
nombre++;
- }
+ }
return std::make_pair(std::string("Echantillon correctement genere : " + lexical_cast<string>(nombre) + " billes"),true);
}
else
=== modified file 'pkg/dem/PreProcessor/SimpleShear.hpp'
--- pkg/dem/PreProcessor/SimpleShear.hpp 2009-12-09 17:11:51 +0000
+++ pkg/dem/PreProcessor/SimpleShear.hpp 2010-02-01 15:20:54 +0000
@@ -63,7 +63,6 @@
void createBox(shared_ptr<Body>& body, Vector3r position, Vector3r extents);
void createSphere(shared_ptr<Body>& body, Vector3r position, Real radius);
void createActors(shared_ptr<Scene>& rootBody);
- void positionRootBody(shared_ptr<Scene>& rootBody);
//method to create a list (containing the positions of centers and radii) of n non interpenetrating spheres, occupying a rectangle with a given (rather high) porosity (issued from TriaxialTest) :
string GenerateCloud(vector<BasicSphere>& sphere_list,Vector3r lowerCorner,Vector3r upperCorner,long number,Real rad_std_dev, Real porosity);
// to create the same list but by reading a text file containing the informations :
Follow ups
