← Back to team overview

yade-dev team mailing list archive

[Branch ~yade-dev/yade/trunk] Rev 2118: I added the (updated) ThreeDTriaxialEngine, and it is also to see if we managed to make bzr worki...

 

------------------------------------------------------------
revno: 2118
committer: Luc Sibille <luc@pc-calc-luc2>
branch nick: trunk
timestamp: Tue 2010-03-30 19:25:38 +0200
message:
  I added the (updated) ThreeDTriaxialEngine, and it is also to see if we managed to make bzr working through an ssh tunnel
added:
  pkg/dem/Engine/PartialEngine/ThreeDTriaxialEngine.cpp
  pkg/dem/Engine/PartialEngine/ThreeDTriaxialEngine.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.
=== added file 'pkg/dem/Engine/PartialEngine/ThreeDTriaxialEngine.cpp'
--- pkg/dem/Engine/PartialEngine/ThreeDTriaxialEngine.cpp	1970-01-01 00:00:00 +0000
+++ pkg/dem/Engine/PartialEngine/ThreeDTriaxialEngine.cpp	2010-03-30 17:25:38 +0000
@@ -0,0 +1,216 @@
+/*************************************************************************
+*  Copyright (C) 2009 by Luc Sibille                                     *
+*  luc.sibille@xxxxxxxxxxxxxx                                            *
+*                                                                        *
+*  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 "ThreeDTriaxialEngine.hpp"
+#include<yade/core/Scene.hpp>
+#include<yade/core/Omega.hpp>
+
+#include<yade/lib-base/Math.hpp>
+#include<boost/lexical_cast.hpp>
+#include<boost/lambda/lambda.hpp>
+#include<yade/pkg-dem/Shop.hpp>
+#include<yade/core/Interaction.hpp>
+#include<yade/pkg-common/Sphere.hpp>
+#include<yade/pkg-dem/FrictPhys.hpp>
+#include<yade/pkg-common/ElastMat.hpp>
+
+class CohesiveFrictionalRelationships;
+
+CREATE_LOGGER(ThreeDTriaxialEngine);
+YADE_PLUGIN((ThreeDTriaxialEngine));
+
+
+ThreeDTriaxialEngine::ThreeDTriaxialEngine()
+{
+	translationAxisy=Vector3r(0,1,0);
+	translationAxisx=Vector3r(1,0,0);
+	translationAxisz=Vector3r(0,0,1);
+	strainRate1=0;
+	currentStrainRate1=0;
+	strainRate2=0;
+	currentStrainRate2=0;
+	strainRate2=0;
+	currentStrainRate2=0;
+	//StabilityCriterion=0.001;
+	UnbalancedForce = 1;
+	Key = "";
+	//Iteration = 0;
+	//testEquilibriumInterval = 20;
+	firstRun=true;
+	frictionAngleDegree = -1;
+	updateFrictionAngle=false;
+
+	stressControl_1=false;
+	stressControl_2=false;
+	stressControl_3=false;
+
+ 	boxVolume=0;
+
+}
+
+ThreeDTriaxialEngine::~ThreeDTriaxialEngine()
+{	
+}
+
+
+void ThreeDTriaxialEngine::action()
+{
+  
+	if ( firstRun )
+	{
+		LOG_INFO ( "First run, will initialize!" );
+
+		if (updateFrictionAngle) setContactProperties(frictionAngleDegree);
+		
+		height0 = height; depth0 = depth; width0 = width;
+
+		if (stressControl_1){
+		wall_right_activated=true; wall_left_activated=true; //are the right walls for direction 1?
+		} else {
+		wall_right_activated=false; wall_left_activated=false;
+		}
+
+		if (stressControl_2){
+		wall_bottom_activated=true; wall_top_activated=true;
+		} else {
+		wall_bottom_activated=false; wall_top_activated=false;
+		}
+
+		if (stressControl_3){
+		wall_front_activated=true; wall_back_activated=true; //are the right walls for direction 3?
+		} else {
+		wall_front_activated=false; wall_back_activated=false;
+		}
+
+		internalCompaction=false;  //is needed to avoid a control for internal compaction by the TriaxialStressController engine
+		
+		//isTriaxialCompression=false;
+		isAxisymetric=false; //is needed to avoid a stress control according the parameter sigma_iso (but according to sigma1, sigma2 and sigma3)
+	
+		firstRun=false;
+
+	}
+
+	
+	/*if ( Omega::instance().getCurrentIteration() % testEquilibriumInterval == 0 )
+	{
+	  //updateParameters ( ncb );
+		computeStressStrain ( ncb );
+
+		UnbalancedForce=ComputeUnbalancedForce ( ncb );
+
+		LOG_INFO("UnbalancedForce="<< UnbalancedForce);
+	
+	}*/
+	
+	
+	Real dt = Omega::instance().getTimeStep();
+
+	if(!stressControl_1)  // control in strain if wanted
+	{
+		if ( currentStrainRate1 != strainRate1 ) currentStrainRate1 += ( strainRate1-currentStrainRate1 ) *0.0003;
+
+		State* p_left=Body::byId(wall_left_id,scene)->state.get();
+		p_left->pos += 0.5*currentStrainRate1*width*translationAxisx*dt;
+		State* p_right=Body::byId(wall_right_id,scene)->state.get();
+		p_right->pos -= 0.5*currentStrainRate1*width*translationAxisx*dt;
+
+	} else {
+
+		if ( currentStrainRate1 != strainRate1 ) currentStrainRate1 += ( strainRate1-currentStrainRate1 ) *0.0003;
+		max_vel1 = 0.5*currentStrainRate1*width;
+	}
+
+	
+	if(!stressControl_2)  // control in strain if wanted
+	{
+		if ( currentStrainRate2 != strainRate2 ) currentStrainRate2 += ( strainRate2-currentStrainRate2 ) *0.0003;
+
+		State* p_bottom=Body::byId(wall_bottom_id,scene)->state.get();
+		p_bottom->pos += 0.5*currentStrainRate2*height*translationAxisy*dt;
+		State* p_top=Body::byId(wall_top_id,scene)->state.get();
+		p_top->pos -= 0.5*currentStrainRate2*height*translationAxisy*dt;
+
+	} else {
+
+		if ( currentStrainRate2 != strainRate2 ) currentStrainRate2 += ( strainRate2-currentStrainRate2 ) *0.0003;
+		max_vel2 = 0.5*currentStrainRate2*height;
+	}
+
+
+	if(!stressControl_3)  // control in strain if wanted
+	{
+		if ( currentStrainRate3 != strainRate3 ) currentStrainRate3 += ( strainRate3-currentStrainRate3 ) *0.0003;
+
+
+		State* p_back=Body::byId(wall_back_id,scene)->state.get();
+		p_back->pos += 0.5*currentStrainRate3*depth*translationAxisz*dt;
+		State* p_front=Body::byId(wall_front_id,scene)->state.get();
+		p_front->pos -= 0.5*currentStrainRate3*depth*translationAxisz*dt;
+
+	} else {
+
+		if ( currentStrainRate3 != strainRate3 ) currentStrainRate3 += ( strainRate3-currentStrainRate3 ) *0.0003;
+		max_vel3 = 0.5*currentStrainRate3*depth;
+	}
+	
+	TriaxialStressController::action(); // this function is called to perform the external stress control
+
+}
+
+void ThreeDTriaxialEngine::setContactProperties(Real frictionDegree)
+{
+	scene = Omega::instance().getScene().get();
+	shared_ptr<BodyContainer>& bodies = scene->bodies;
+	FOREACH(const shared_ptr<Body>& b,*scene->bodies){
+		if (b->isDynamic)
+		YADE_PTR_CAST<FrictMat> (b->material)->frictionAngle = frictionDegree * Mathr::PI/180.0;
+	}
+			
+// 	BodyContainer::iterator bi = bodies->begin();
+// 	BodyContainer::iterator biEnd = bodies->end();
+// 	
+// 	for ( ; bi!=biEnd; ++bi)	
+// 	{	
+// 		shared_ptr<Body> b = *bi;
+// 		if (b->isDynamic)
+// 
+// 		YADE_PTR_CAST<FrictMat> (b->material)->frictionAngle = frictionDegree * Mathr::PI/180.0;
+// 	}
+	
+	FOREACH(const shared_ptr<Interaction>& ii, *scene->interactions){
+		if (!ii->isReal()) continue;
+		const shared_ptr<FrictMat>& sdec1 = YADE_PTR_CAST<FrictMat>((*bodies)[(body_id_t) ((ii)->getId1())]->material);
+		const shared_ptr<FrictMat>& sdec2 = YADE_PTR_CAST<FrictMat>((*bodies)[(body_id_t) ((ii)->getId2())]->material);
+		//FIXME - why dynamic_cast fails here?
+		FrictPhys* contactPhysics = YADE_CAST<FrictPhys*>((ii)->interactionPhysics.get());
+		Real fa = sdec1->frictionAngle;
+		Real fb = sdec2->frictionAngle;
+		contactPhysics->frictionAngle			= std::min(fa,fb);
+		contactPhysics->tangensOfFrictionAngle		= std::tan(contactPhysics->frictionAngle);
+	}
+  
+// 	InteractionContainer::iterator ii    = ncb->interactions->begin();
+// 	InteractionContainer::iterator iiEnd = ncb->interactions->end(); 
+// 	for(  ; ii!=iiEnd ; ++ii ) {
+// 		if (!(*ii)->isReal()) continue;
+// 		
+// 		const shared_ptr<FrictMat>& sdec1 = YADE_PTR_CAST<FrictMat>((*bodies)[(body_id_t) ((*ii)->getId1())]->material);
+// 		const shared_ptr<FrictMat>& sdec2 = YADE_PTR_CAST<FrictMat>((*bodies)[(body_id_t) ((*ii)->getId2())]->material);
+// 		
+// 		const shared_ptr<FrictPhys>& contactPhysics = static_pointer_cast<FrictPhys>((*ii)->interactionPhysics);
+// 
+// 		Real fa 	= sdec1->frictionAngle;
+// 		Real fb 	= sdec2->frictionAngle;
+// 
+// 		contactPhysics->frictionAngle			= std::min(fa,fb);
+// 		contactPhysics->tangensOfFrictionAngle		= std::tan(contactPhysics->frictionAngle); 
+//	}
+} 
+
+

=== added file 'pkg/dem/Engine/PartialEngine/ThreeDTriaxialEngine.hpp'
--- pkg/dem/Engine/PartialEngine/ThreeDTriaxialEngine.hpp	1970-01-01 00:00:00 +0000
+++ pkg/dem/Engine/PartialEngine/ThreeDTriaxialEngine.hpp	2010-03-30 17:25:38 +0000
@@ -0,0 +1,88 @@
+/*************************************************************************
+*  Copyright (C) 2009 by Luc Sibille                                     *
+*  luc.sibille@xxxxxxxxxxxxxx                                            *
+*                                                                        *
+*  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/PartialEngine.hpp>
+#include<yade/lib-base/Math.hpp>
+#include<yade/pkg-dem/TriaxialStressController.hpp>
+#include<string>
+
+
+
+/** \brief Class for controlling in stress or in strain with respect to each spatial direction a cubical assembly of particles.
+ *
+ * The engine perform a triaxial compression with a control in direction "i" in stress "if (stressControl_i)" else in strain.
+ * For a stress control the imposed stress is specified by "sigma_i" with a "max_veli" depending on "strainRatei". To obtain the same strain rate in stress control than in strain control you need to set "wallDamping = 0.8".
+ * For a strain control the imposed strain is specified by "strainRatei".
+ *
+ */
+
+class ThreeDTriaxialEngine : public TriaxialStressController
+{
+	public :
+		ThreeDTriaxialEngine();
+		virtual ~ThreeDTriaxialEngine();
+				
+		//! Strain velocity (./s)
+		Real strainRate1;
+		Real currentStrainRate1;
+		Real strainRate2;
+		Real currentStrainRate2;
+		Real strainRate3;
+		Real currentStrainRate3;
+		//! Max ratio of resultant forces on mean contact force
+		Real UnbalancedForce;
+		//! Value of UnbalancedForce for which the system is considered stable
+		//Real StabilityCriterion;
+		//! Value of axial deformation for which the simulation must stop
+		//Real epsilonMax;
+		//! Current value of axial deformation during confined loading (is reference to strain[1])
+		//Real& uniaxialEpsilonCurr;
+		//! Value of friction to use in the test "if (updateFrictionAngle)"
+		Real frictionAngleDegree;
+		//! Swith to activate if an update of the intergranular friction angle is required
+		bool updateFrictionAngle;
+
+		//! switches to choose a stress or a strain control in directions 1, 2 or 3
+		bool stressControl_1;
+		bool stressControl_2;
+		bool stressControl_3;
+		
+
+		Vector3r translationAxisy;
+		Vector3r translationAxisx;
+		Vector3r translationAxisz;
+
+		//! is this the beginning of the simulation, after reading the scene? -> it is the first time that Yade passes trought the engine ThreeDTriaxialEngine
+		bool firstRun;
+		//int FinalIterationPhase1, Iteration, testEquilibriumInterval;
+		//int Iteration, testEquilibriumInterval;
+		
+		std::string Key;//A code that is appended to file names to help distinguish between different simulations
+				
+		virtual void action();
+
+		
+		///Change physical properties of interactions and/or bodies in the middle of a simulation (change only friction for the moment, complete this function to set cohesion and others before compression test)
+		void setContactProperties(Real frictionDegree);
+
+		
+
+	protected :
+		REGISTER_ATTRIBUTES(TriaxialStressController,(strainRate1)(currentStrainRate1)(strainRate2)(currentStrainRate2)(strainRate3)(currentStrainRate3)/*(Phase1)*/(UnbalancedForce)/*(StabilityCriterion)*//*(translationAxis)(compressionActivated)(autoCompressionActivation)(autoStopSimulation)*//*(testEquilibriumInterval)*//*(currentState)(previousState)(sigmaIsoCompaction)(previousSigmaIso)(sigmaLateralConfinement)*/(Key)(frictionAngleDegree)(updateFrictionAngle)	/*(epsilonMax)*/ /*(uniaxialEpsilonCurr)*//*(isotropicCompaction)*/(spheresVolume)/*(fixedPorosity)*/(stressControl_1)(stressControl_2)(stressControl_3)(sigma1)(sigma2)(sigma3));
+		
+	REGISTER_CLASS_NAME(ThreeDTriaxialEngine);
+	REGISTER_BASE_CLASS_NAME(TriaxialStressController);
+	DECLARE_LOGGER;
+};
+
+REGISTER_SERIALIZABLE(ThreeDTriaxialEngine);
+
+


Follow ups