← Back to team overview

yade-dev team mailing list archive

RE : Modified GlobalStiffnessTimeStepper for visco-elastic contact law

 

Hi all, 
Here is the final version of GlobalStiffnessTimeStepper that I modified to take into account the viscous contribution when using Law2_ScGeom_ViscElPhys_Basic. 
I didn't modify the elastic part, so that it will work exactly the same for people using it as before (considering only elastic contribution). I just added an option : viscEl which, when True, allows to calculate the time step for ViscElPhys as described before.
I tested it and everything works well.

I put the patch below for GlobalStiffnessTimeStepper.cpp and GlobalStiffnessTimeStepper.hpp, can someone add the change on github ? 
Thanks

Raphael



GlobalStiffnessTimeStepper.cpp



diff --git a/pkg/dem/GlobalStiffnessTimeStepper.cpp b/pkg/dem/GlobalStiffnessTimeStepper.cpp
index 8a41ceb..0205655 100644
--- a/pkg/dem/GlobalStiffnessTimeStepper.cpp
+++ b/pkg/dem/GlobalStiffnessTimeStepper.cpp
@@ -14,6 +14,7 @@
 #include<yade/core/Scene.hpp>
 #include<yade/core/Clump.hpp>
 #include<yade/pkg/dem/Shop.hpp>
+#include<yade/pkg/dem/ViscoelasticPM.hpp>
 
 CREATE_LOGGER(GlobalStiffnessTimeStepper);
 YADE_PLUGIN((GlobalStiffnessTimeStepper));
@@ -25,13 +26,16 @@ void GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& bo
 	State* sdec=body->state.get();
 	Vector3r&  stiffness= stiffnesses[body->getId()];
 	Vector3r& Rstiffness=Rstiffnesses[body->getId()];
-	
 	if(body->isClump()) {// if clump, we sum stifnesses of all members
 		const shared_ptr<Clump>& clump=YADE_PTR_CAST<Clump>(body->shape);
 		FOREACH(Clump::MemberMap::value_type& B, clump->members){
 			const shared_ptr<Body>& b = Body::byId(B.first,scene);
 			stiffness+=stiffnesses[b->getId()];
 			Rstiffness+=Rstiffnesses[b->getId()];
+			if (viscEl == true){
+				viscosities[body->getId()]+=viscosities[b->getId()];
+				Rviscosities[body->getId()]+=Rviscosities[b->getId()];
+			}
 		}
 	}
 	
@@ -39,13 +43,15 @@ void GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& bo
 		if (densityScaling) sdec->densityScaling = min(1.0001*sdec->densityScaling, timestepSafetyCoefficient*pow(defaultDt/targetDt,2.0));
 		return; // not possible to compute!
 	}
-	
+
+	//Determine the elastic minimum eigenperiod (and if required determine also the viscous one separately and take the minimum of the two)
+
+	//Elastic
 	Real dtx, dty, dtz;
 	Real dt = max( max (stiffness.x(), stiffness.y()), stiffness.z() );
 	if (dt!=0) {
 		dt = sdec->mass/dt;  computedSomething = true;}//dt = squared eigenperiod of translational motion 
 	else dt = Mathr::MAX_REAL;
-	
 	if (Rstiffness.x()!=0) {
 		dtx = sdec->inertia.x()/Rstiffness.x();  computedSomething = true;}//dtx = squared eigenperiod of rotational motion around x
 	else dtx = Mathr::MAX_REAL;	
@@ -56,15 +62,43 @@ void GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& bo
 		dtz = sdec->inertia.z()/Rstiffness.z();  computedSomething = true;}
 	else dtz = Mathr::MAX_REAL;
 	
-	Real Rdt =  std::min( std::min (dtx, dty), dtz );//Rdt = smallest squared eigenperiod for rotational motions
+	Real Rdt =  std::min( std::min (dtx, dty), dtz );//Rdt = smallest squared eigenperiod for elastic rotational motions
 	dt = 1.41044*timestepSafetyCoefficient*std::sqrt(std::min(dt,Rdt));//1.41044 = sqrt(2)
+	
+	//Viscous 
+	if (viscEl == true){		
+		Vector3r&  viscosity = viscosities[body->getId()];
+		Vector3r& Rviscosity = Rviscosities[body->getId()];
+		Real dtx_visc, dty_visc, dtz_visc;
+		Real dt_visc = max(max(viscosity.x(), viscosity.y()), viscosity.z() );
+		if (dt_visc!=0) {
+			dt_visc = sdec->mass/dt_visc;  computedSomething = true;}//dt = eigenperiod of the viscous translational motion 
+		else {dt_visc = Mathr::MAX_REAL;}
+
+		if (Rviscosity.x()!=0) {
+			dtx_visc = sdec->inertia.x()/Rviscosity.x();  computedSomething = true;}//dtx = eigenperiod of viscous rotational motion around x
+		else dtx_visc = Mathr::MAX_REAL;	
+		if (Rviscosity.y()!=0) {
+			dty_visc = sdec->inertia.y()/Rviscosity.y();  computedSomething = true;}
+		else dty_visc = Mathr::MAX_REAL;
+		if (Rviscosity.z()!=0) {
+			dtz_visc = sdec->inertia.z()/Rviscosity.z();  computedSomething = true;}
+		else dtz_visc = Mathr::MAX_REAL;
+	
+		Real Rdt_visc =  std::min( std::min (dtx_visc, dty_visc), dtz_visc );//Rdt = smallest squared eigenperiod for viscous rotational motions
+		dt_visc = 2*timestepSafetyCoefficient*std::min(dt_visc,Rdt_visc);
+
+		//Take the minimum between the elastic and viscous minimum eigenperiod. 
+		dt = std::min(dt,dt_visc);
+	}
+
 	//if there is a target dt, then we apply density scaling on the body, the inertia used in Newton will be mass*scaling, the weight is unmodified
 	if (densityScaling) {
 		sdec->densityScaling = min(sdec->densityScaling,timestepSafetyCoefficient*pow(dt /targetDt,2.0));
 		newDt=targetDt;
 	}
 	//else we update dt normaly
-	else {newDt = std::min(dt,newDt);}
+	else {newDt = std::min(dt,newDt);}   
 }
 
 bool GlobalStiffnessTimeStepper::isActivated()
@@ -97,7 +131,7 @@ void GlobalStiffnessTimeStepper::computeTimeStep(Scene* ncb)
 		computedOnce = true;}
 	else if (!computedOnce) scene->dt=defaultDt;
 	LOG_INFO("computed timestep " << newDt <<
-			(scene->dt==newDt ? string(", appplied") :
+			(scene->dt==newDt ? string(", applied") :
 			string(", BUT timestep is ")+lexical_cast<string>(scene->dt))<<".");
 }
 
@@ -107,15 +141,24 @@ void GlobalStiffnessTimeStepper::computeStiffnesses(Scene* rb){
 	if(size<rb->bodies->size()){
 		size=rb->bodies->size();
 		stiffnesses.resize(size); Rstiffnesses.resize(size);
+		if (viscEl == true){
+			viscosities.resize(size); Rviscosities.resize(size);
+			}
 	}
 	/* reset stored values */
 	memset(& stiffnesses[0],0,sizeof(Vector3r)*size);
 	memset(&Rstiffnesses[0],0,sizeof(Vector3r)*size);
+	if (viscEl == true){
+		memset(& viscosities[0],0,sizeof(Vector3r)*size);
+		memset(&Rviscosities[0],0,sizeof(Vector3r)*size);
+	}
+
 	FOREACH(const shared_ptr<Interaction>& contact, *rb->interactions){
 		if(!contact->isReal()) continue;
 
 		GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(contact->geom.get()); assert(geom);
 		NormShearPhys* phys=YADE_CAST<NormShearPhys*>(contact->phys.get()); assert(phys);
+
 		// all we need for getting stiffness
 		Vector3r& normal=geom->normal; Real& kn=phys->kn; Real& ks=phys->ks; Real& radius1=geom->refR1; Real& radius2=geom->refR2;
 		Real fn = (static_cast<NormShearPhys *> (contact->phys.get()))->normalForce.squaredNorm();
@@ -134,10 +177,35 @@ void GlobalStiffnessTimeStepper::computeStiffnesses(Scene* rb){
 				std::pow(normal.x(),2)+std::pow(normal.z(),2),
 				std::pow(normal.x(),2)+std::pow(normal.y(),2));
 		diag_Rstiffness *= ks;
+
+
 		//NOTE : contact laws with moments would be handled correctly by summing directly bending+twisting stiffness to diag_Rstiffness. The fact that there is no problem currently with e.g. cohesiveFrict law is probably because final computed dt is constrained by translational motion, not rotations.
 		stiffnesses [contact->getId1()]+=diag_stiffness;
 		Rstiffnesses[contact->getId1()]+=diag_Rstiffness*pow(radius1,2);
 		stiffnesses [contact->getId2()]+=diag_stiffness;
-		Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);
+		Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);	
+
+		//Same for the Viscous part, if required
+		if (viscEl == true){
+			ViscElPhys* viscPhys = YADE_CAST<ViscElPhys*>(contact->phys.get()); assert(viscPhys);
+			Real& cn=viscPhys->cn; Real& cs=viscPhys->cs;
+			//Diagonal terms of the translational viscous matrix
+			Vector3r diag_viscosity = Vector3r(std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));
+			diag_viscosity *= cn-cs;
+			diag_viscosity = diag_viscosity + Vector3r(1,1,1)*cs;
+			//diagonal terms of the rotational viscous matrix
+			Vector3r diag_Rviscosity =
+				Vector3r(std::pow(normal.y(),2)+std::pow(normal.z(),2),
+					std::pow(normal.x(),2)+std::pow(normal.z(),2),
+					std::pow(normal.x(),2)+std::pow(normal.y(),2));
+			diag_Rviscosity *= cs;			
+			
+			// Add the contact stiffness matrix to the two particles one
+			viscosities [contact->getId1()]+=diag_viscosity;
+			Rviscosities[contact->getId1()]+=diag_Rviscosity*pow(radius1,2);
+			viscosities [contact->getId2()]+=diag_viscosity;
+			Rviscosities[contact->getId2()]+=diag_Rviscosity*pow(radius2,2);
+		}
+
 	}
 }






GlobalStiffnessTimeStepper.hpp






diff --git a/pkg/dem/GlobalStiffnessTimeStepper.hpp b/pkg/dem/GlobalStiffnessTimeStepper.hpp
index 1d60f88..bd657fa 100644
--- a/pkg/dem/GlobalStiffnessTimeStepper.hpp
+++ b/pkg/dem/GlobalStiffnessTimeStepper.hpp
@@ -24,6 +24,8 @@ class GlobalStiffnessTimeStepper : public TimeStepper
 	private :
 		vector<Vector3r> stiffnesses;
 		vector<Vector3r> Rstiffnesses;
+		vector<Vector3r> viscosities;
+		vector<Vector3r> Rviscosities;
 		void computeStiffnesses(Scene*);
 
 		Real		newDt;
@@ -37,13 +39,14 @@ class GlobalStiffnessTimeStepper : public TimeStepper
 		virtual void computeTimeStep(Scene*);
 		virtual bool isActivated();
 		YADE_CLASS_BASE_DOC_ATTRS_CTOR(
-			GlobalStiffnessTimeStepper,TimeStepper,"An engine assigning the time-step as a fraction of the minimum eigen-period in the problem. The derivation is detailed in the chapter on `DEM formulation <https://www.yade-dem.org/doc/formulation.html#dem-simulations>`_",
+			GlobalStiffnessTimeStepper,TimeStepper,"An engine assigning the time-step as a fraction of the minimum eigen-period in the problem. The derivation is detailed in the chapter on `DEM formulation <https://www.yade-dem.org/doc/formulation.html#dem-simulations>`_. The viscEl option enables to evaluate the timestep in a similar way for the visco-elastic contact law :yref:`Law2_ScGeom_ViscElPhys_Basic`, more detail in :yref:`GlobalStiffnessTimestepper::viscEl`. ",
 			((Real,defaultDt,-1,,"used as the initial value of the timestep (especially useful in the first steps when no contact exist). If negative, it will be defined by :yref:`utils.PWaveTimeStep` * :yref:`GlobalStiffnessTimeStepper::timestepSafetyCoefficient`"))
 			((Real,maxDt,Mathr::MAX_REAL,,"if positive, used as max value of the timestep whatever the computed value"))
 			((Real,previousDt,1,,"last computed dt |yupdate|"))
 			((Real,timestepSafetyCoefficient,0.8,,"safety factor between the minimum eigen-period and the final assigned dt (less than 1))"))
 			((bool,densityScaling,false,,"|yupdate| don't modify this value if you don't plan to modify the scaling factor manually for some bodies. In most cases, it is enough to set :yref:`NewtonIntegrator::densityScaling` and let this one be adjusted automatically."))
-			((Real,targetDt,1,,"if :yref:`NewtonIntegrator::densityScaling` is active, this value will be used as the simulation  timestep and the scaling will use this value of dt as the target value. The value of targetDt is arbitrary and should have no effect in the result in general. However if some bodies have imposed velocities, for instance, they will move more or less per each step depending on this value.")),
+			((Real,targetDt,1,,"if :yref:`NewtonIntegrator::densityScaling` is active, this value will be used as the simulation  timestep and the scaling will use this value of dt as the target value. The value of targetDt is arbitrary and should have no effect in the result in general. However if some bodies have imposed velocities, for instance, they will move more or less per each step depending on this value."))
+			((bool,viscEl,false,,"To use with :yref:`ViscElPhys`. if True, evaluate separetly the minimum eigen-period in the problem considering only the elastic contribution on one hand (spring only), and only the viscous contribution on the other hand (dashpot only). Take then the minimum of the two and use the safety coefficient :yref:`GlobalStiffnessTimestepper::timestepSafetyCoefficient` to take into account the possible coupling between the two contribution.")),
 			computedOnce=false;)
 		DECLARE_LOGGER;
 };










-------- Message d'origine--------
De: Raphaël Maurin
Date: mar. 03/09/2013 12:21
À: yade-dev@xxxxxxxxxxxxxxxxxxx
Objet : Modified GlobalStiffnessTimeStepper for visco-elastic contact law
 
Hi all, 

I developed a "copy" of GlobalStiffnessTimeStepper which evaluate the time step of a system of particle with Law2_ScGeom_ViscElPhys_Basic contact law (spring + viscous damping).

The principle of the time step evaluation is : 
We consider that the elastic and the viscous degrees of freedom are decoupled. From that, we evaluate the critical time step of the system of particles considering the elastic part only (dt_el) on one hand, and considering the viscous part only (dt_visc) on the other hand. Then, we take the minimum of the two timestep and multiply it by a safety factor to account for the possible coupling between the two effect (dt = alpha*min(dt_el,dt_visc). 

To do so, I took GlobalStiffnessTimeStepper (GSTS) which is already evaluating the elastic critical time step, and implemented similarly a calculation for the viscous critical time step of the system from the the viscosity matrix of each particles. The changes are not so important between the version I implemented and the last one (see the code below). I tested it in cases where the normal viscous part is dominant and the new version gives a lower time step which allows to keep the caculus stable, while in the old version everything explode, so it seems to work.

I have some questions about that :
- Does the evaluation of the time step seems reasonable to you? 
- What can I do to test it quantitavely ? (if I try a collision between two beads, GSTS does not work..)
- To add it to the code, should I add it as a new "file"/"function" or should I implement it as an option of GSTS ? 

And about GSTS in general : 
- why is GSTS not working when we consider a gravity deposition or a collision ?? (even when defaultDt is small enough)

Thanks, 

Raphael



GlobalStiffnessTimeStepper.cpp

/*************************************************************************
*  Copyright (C) 2006 by Bruno Chareyre                                  *
*  bruno.chareyre@xxxxxxxxxxx                                            *
*                                                                        *
*  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"GlobalStiffnessTimeStepper.hpp"
#include<yade/pkg/dem/FrictPhys.hpp>
#include<yade/pkg/dem/ScGeom.hpp>
#include<yade/pkg/dem/DemXDofGeom.hpp>
#include<yade/core/Interaction.hpp>
#include<yade/core/Scene.hpp>
#include<yade/core/Clump.hpp>
#include<yade/pkg/dem/Shop.hpp>
#include<yade/pkg/dem/ViscoelasticPM.hpp>

CREATE_LOGGER(GlobalStiffnessTimeStepper);
YADE_PLUGIN((GlobalStiffnessTimeStepper));

GlobalStiffnessTimeStepper::~GlobalStiffnessTimeStepper() {}

void GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& body, Scene * ncb)
{
	State* sdec=body->state.get();
	Vector3r&  stiffness= stiffnesses[body->getId()];
	Vector3r& Rstiffness=Rstiffnesses[body->getId()];
	Vector3r& viscosity= viscosities[body->getId()];
	Vector3r& Rviscosity= Rviscosities[body->getId()];

	if(body->isClump()) {// if clump, we sum stifnesses of all members
		const shared_ptr<Clump>& clump=YADE_PTR_CAST<Clump>(body->shape);
		FOREACH(Clump::MemberMap::value_type& B, clump->members){
			const shared_ptr<Body>& b = Body::byId(B.first,scene);
			stiffness+=stiffnesses[b->getId()];
			Rstiffness+=Rstiffnesses[b->getId()];
			viscosity+=viscosities[b->getId()];
			Rviscosity+=Rviscosities[b->getId()];
		}
	}
	
	if(!sdec || stiffness==Vector3r::Zero()){
		if (densityScaling) sdec->densityScaling = min(1.0001*sdec->densityScaling, timestepSafetyCoefficient*pow(defaultDt/targetDt,2.0));
		return; // not possible to compute!
	}

	//Determine the elastic critical time step and the viscous one separately (as if they were decoupled) and take the minimum	time step with a safety coefficient. 

	//Elastic
	Real dtx, dty, dtz;
	Real dt = max( max (stiffness.x(), stiffness.y()), stiffness.z() );
	if (dt!=0) {
		dt = sdec->mass/dt;  computedSomething = true;}//dt = squared eigenperiod of translational motion 
	else dt = Mathr::MAX_REAL;
	if (Rstiffness.x()!=0) {
		dtx = sdec->inertia.x()/Rstiffness.x();  computedSomething = true;}//dtx = squared eigenperiod of rotational motion around x
	else dtx = Mathr::MAX_REAL;	
	if (Rstiffness.y()!=0) {
		dty = sdec->inertia.y()/Rstiffness.y();  computedSomething = true;}
	else dty = Mathr::MAX_REAL;
	if (Rstiffness.z()!=0) {
		dtz = sdec->inertia.z()/Rstiffness.z();  computedSomething = true;}
	else dtz = Mathr::MAX_REAL;
	
	Real Rdt =  std::min( std::min (dtx, dty), dtz );//Rdt = smallest squared eigenperiod for elastic rotational motions
	dt = 1.41044*timestepSafetyCoefficient*std::sqrt(std::min(dt,Rdt));//1.41044 = sqrt(2)
	
	//Viscous
	Real dtx_visc, dty_visc, dtz_visc;
	Real dt_visc = max(max(viscosity.x(), viscosity.y()), viscosity.z() );
	if (dt_visc!=0) {
		dt_visc = sdec->mass/dt_visc;  computedSomething = true;}//dt = eigenperiod of the viscous translational motion 
	else {dt_visc = Mathr::MAX_REAL;}

	if (Rviscosity.x()!=0) {
		dtx_visc = sdec->inertia.x()/Rviscosity.x();  computedSomething = true;}//dtx = eigenperiod of viscous rotational motion around x
	else dtx_visc = Mathr::MAX_REAL;	
	if (Rviscosity.y()!=0) {
		dty_visc = sdec->inertia.y()/Rviscosity.y();  computedSomething = true;}
	else dty_visc = Mathr::MAX_REAL;
	if (Rviscosity.z()!=0) {
		dtz_visc = sdec->inertia.z()/Rviscosity.z();  computedSomething = true;}
	else dtz_visc = Mathr::MAX_REAL;
	
	Real Rdt_visc =  std::min( std::min (dtx_visc, dty_visc), dtz_visc );//Rdt = smallest squared eigenperiod for viscous rotational motions
	dt_visc = 2*timestepSafetyCoefficient*std::min(dt_visc,Rdt_visc);

	//Take the minimum between the elastic and viscous critical time step. 
	dt = std::min(dt,dt_visc);

	//if there is a target dt, then we apply density scaling on the body, the inertia used in Newton will be mass*scaling, the weight is unmodified
	if (densityScaling) {
		sdec->densityScaling = min(sdec->densityScaling,timestepSafetyCoefficient*pow(dt /targetDt,2.0));
		newDt=targetDt;
	}
	//else we update dt normaly
	else {newDt = std::min(dt,newDt);}   
}

bool GlobalStiffnessTimeStepper::isActivated()
{
	return (active && ((!computedOnce) || (scene->iter % timeStepUpdateInterval == 0) || (scene->iter < (long int) 2) ));
}

void GlobalStiffnessTimeStepper::computeTimeStep(Scene* ncb)
{
	// for some reason, this line is necessary to have correct functioning (no idea _why_)
	// see scripts/test/compare-identical.py, run with or without active=active.
	active=active;
	if (defaultDt<0) defaultDt= timestepSafetyCoefficient*Shop::PWaveTimeStep(Omega::instance().getScene());
	computeStiffnesses(ncb);

	shared_ptr<BodyContainer>& bodies = ncb->bodies;
	newDt = Mathr::MAX_REAL;
	computedSomething=false;
	BodyContainer::iterator bi    = bodies->begin();
	BodyContainer::iterator biEnd = bodies->end();
	for(  ; bi!=biEnd ; ++bi ){
		shared_ptr<Body> b = *bi;
		if (b->isDynamic() && !b->isClumpMember()) findTimeStepFromBody(b, ncb);
		
	}
	if(densityScaling) (newDt=targetDt);
	if(computedSomething || densityScaling){
		previousDt = min ( min(newDt , maxDt), 1.05*previousDt );// at maximum, dt will be multiplied by 1.05 in one iterration, this is to prevent brutal switches from 0.000... to 1 in some computations
		scene->dt=previousDt;
		computedOnce = true;}
	else if (!computedOnce) scene->dt=defaultDt;
	LOG_INFO("computed timestep " << newDt <<
			(scene->dt==newDt ? string(", applied") :
			string(", BUT timestep is ")+lexical_cast<string>(scene->dt))<<".");
}

void GlobalStiffnessTimeStepper::computeStiffnesses(Scene* rb){
	/* check size */
	size_t size=stiffnesses.size();
	if(size<rb->bodies->size()){
		size=rb->bodies->size();
		stiffnesses.resize(size); Rstiffnesses.resize(size);
		viscosities.resize(size); Rviscosities.resize(size);
	}
	/* reset stored values */
	memset(& stiffnesses[0],0,sizeof(Vector3r)*size);
	memset(&Rstiffnesses[0],0,sizeof(Vector3r)*size);
	memset(& viscosities[0],0,sizeof(Vector3r)*size);
	memset(& Rviscosities[0],0,sizeof(Vector3r)*size);

	FOREACH(const shared_ptr<Interaction>& contact, *rb->interactions){
		if(!contact->isReal()) continue;

		GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(contact->geom.get()); assert(geom);
		ViscElPhys* phys = YADE_CAST<ViscElPhys*>(contact->phys.get()); assert(phys);
		// all we need for getting stiffness
		Vector3r& normal=geom->normal; Real& kn=phys->kn; Real& ks=phys->ks; Real& radius1=geom->refR1; Real& radius2=geom->refR2;
		Real& cn=phys->cn; Real& cs=phys->cs;
		Real fn = (static_cast<NormShearPhys *> (contact->phys.get()))->normalForce.squaredNorm();
		if (fn==0) continue;//Is it a problem with some laws? I still don't see why.
		
		//Diagonal terms of the translational stiffness matrix
		Vector3r diag_stiffness = Vector3r(std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));
		diag_stiffness *= kn-ks;
		diag_stiffness = diag_stiffness + Vector3r(1,1,1)*ks;

		//diagonal terms of the rotational stiffness matrix
		// Vector3r branch1 = currentContactGeometry->normal*currentContactGeometry->radius1;
		// Vector3r branch2 = currentContactGeometry->normal*currentContactGeometry->radius2;
		Vector3r diag_Rstiffness =
			Vector3r(std::pow(normal.y(),2)+std::pow(normal.z(),2),
				std::pow(normal.x(),2)+std::pow(normal.z(),2),
				std::pow(normal.x(),2)+std::pow(normal.y(),2));
		diag_Rstiffness *= ks;

		//Diagonal terms of the translational viscous matrix
		Vector3r diag_viscosity = Vector3r(std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));
		diag_viscosity *= cn-cs;
		diag_viscosity = diag_viscosity + Vector3r(1,1,1)*cs;
		//diagonal terms of the rotational viscous matrix
		Vector3r diag_Rviscosity =
			Vector3r(std::pow(normal.y(),2)+std::pow(normal.z(),2),
				std::pow(normal.x(),2)+std::pow(normal.z(),2),
				std::pow(normal.x(),2)+std::pow(normal.y(),2));
		diag_Rviscosity *= cs;

		//NOTE : contact laws with moments would be handled correctly by summing directly bending+twisting stiffness to diag_Rstiffness. The fact that there is no problem currently with e.g. cohesiveFrict law is probably because final computed dt is constrained by translational motion, not rotations.
		stiffnesses [contact->getId1()]+=diag_stiffness;
		Rstiffnesses[contact->getId1()]+=diag_Rstiffness*pow(radius1,2);
		viscosities [contact->getId1()]+=diag_viscosity;
		Rviscosities[contact->getId1()]+=diag_Rviscosity*pow(radius1,2);
		stiffnesses [contact->getId2()]+=diag_stiffness;
		Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);
		viscosities [contact->getId2()]+=diag_viscosity;
		Rviscosities[contact->getId2()]+=diag_Rviscosity*pow(radius2,2);
	}
}











GlobalStiffnessTimeStepper.hpp



/*************************************************************************
*  Copyright (C) 2006 by Bruno Chareyre                                  *
*  bruno.chareyre@xxxxxxxxxxx                                            *
*                                                                        *
*  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/TimeStepper.hpp>


/*! \brief Compute the critical timestep of the leap-frog scheme based on global stiffness of bodies
	See usage details in TriaxialTest documentation (TriaxialTest is also a good example of how to use this class)
 */

class Interaction;
class BodyContainer;
class Scene;

class GlobalStiffnessTimeStepper : public TimeStepper
{
	private :
		vector<Vector3r> stiffnesses;
		vector<Vector3r> Rstiffnesses;
		vector<Vector3r> viscosities;
		vector<Vector3r> Rviscosities;
		void computeStiffnesses(Scene*);

		Real		newDt;
		bool		computedSomething,
 				computedOnce;
		void findTimeStepFromBody(const shared_ptr<Body>& body, Scene * ncb);
	
	public :
		virtual ~GlobalStiffnessTimeStepper();
	
		virtual void computeTimeStep(Scene*);
		virtual bool isActivated();
		YADE_CLASS_BASE_DOC_ATTRS_CTOR(
			GlobalStiffnessTimeStepper,TimeStepper,"An engine assigning the time-step as a fraction of the minimum eigen-period in the problem. The derivation is detailed in the chapter on `DEM formulation <https://www.yade-dem.org/doc/formulation.html#dem-simulations>`_",
			((Real,defaultDt,-1,,"used as the initial value of the timestep (especially useful in the first steps when no contact exist). If negative, it will be defined by :yref:`utils.PWaveTimeStep` * :yref:`GlobalStiffnessTimeStepper::timestepSafetyCoefficient`"))
			((Real,maxDt,Mathr::MAX_REAL,,"if positive, used as max value of the timestep whatever the computed value"))
			((Real,previousDt,1,,"last computed dt |yupdate|"))
			((Real,timestepSafetyCoefficient,0.8,,"safety factor between the minimum eigen-period and the final assigned dt (less than 1))"))
			((bool,densityScaling,false,,"|yupdate| don't modify this value if you don't plan to modify the scaling factor manually for some bodies. In most cases, it is enough to set :yref:`NewtonIntegrator::densityScaling` and let this one be adjusted automatically."))
			((Real,targetDt,1,,"if :yref:`NewtonIntegrator::densityScaling` is active, this value will be used as the simulation  timestep and the scaling will use this value of dt as the target value. The value of targetDt is arbitrary and should have no effect in the result in general. However if some bodies have imposed velocities, for instance, they will move more or less per each step depending on this value.")),
			computedOnce=false;)
		DECLARE_LOGGER;
};

REGISTER_SERIALIZABLE(GlobalStiffnessTimeStepper);





References