yade-dev team mailing list archive

[Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>]

>  If you agree could you commit it, please? Thanks.
Oh, sorry, I forgot that point... :-/
I'm sending you the sources, because I have also diffs for the OPenMP 
plastic accumulator and I need to compile before commiting (I'm already 
compiling another trunk and it will take time).
I will commit almost the same files tomorrow if they compile (they 
should), but OpenMP thing will be a bit more clean.

Bruno

/*************************************************************************
*  Copyright (C) 2004 by Olivier Galizzi                                 *
*  olivier.galizzi@xxxxxxx                                               *
*                                                                        *
*  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"ElasticContactLaw.hpp"
#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/pkg-dem/FrictPhys.hpp>
#include<yade/pkg-dem/DemXDofGeom.hpp>
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
#include<yade/core/Scene.hpp>

YADE_PLUGIN((Law2_ScGeom_FrictPhys_Basic)(Law2_Dem3DofGeom_FrictPhys_Basic)(ElasticContactLaw)(Law2_Dem6DofGeom_FrictPhys_Beam));

Real Law2_ScGeom_FrictPhys_Basic::Real0=0;

void ElasticContactLaw::action()
{
	if(!functor) functor=shared_ptr<Law2_ScGeom_FrictPhys_Basic>(new Law2_ScGeom_FrictPhys_Basic);
	functor->sdecGroupMask=sdecGroupMask;
	functor->useShear=useShear;
	functor->neverErase=neverErase;
	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
		if(!I->isReal()) continue;
		#ifdef YADE_DEBUG
			// these checks would be redundant in the functor (LawDispatcher does that already)
			if(!dynamic_cast<ScGeom*>(I->interactionGeometry.get()) || !dynamic_cast<FrictPhys*>(I->interactionPhysics.get())) continue;	
		#endif
			functor->go(I->interactionGeometry, I->interactionPhysics, I.get(), scene);
	}
}

Real Law2_ScGeom_FrictPhys_Basic::elasticEnergy()
{
	Real energy=0;
	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
		if(!I->isReal()) continue;
		FrictPhys* phys = dynamic_cast<FrictPhys*>(I->interactionPhysics.get());
		if(phys) {
			energy += 0.5*(phys->normalForce.squaredNorm()/phys->kn + phys->shearForce.squaredNorm()/phys->ks);}
	}
	return energy;
}

Real Law2_ScGeom_FrictPhys_Basic::getPlasticDissipation() {return (Real) plasticDissipation;}

void Law2_ScGeom_FrictPhys_Basic::initPlasticDissipation(Real initVal) {plasticDissipation.reset(); plasticDissipation+=initVal;}

CREATE_LOGGER(Law2_ScGeom_FrictPhys_Basic);
void Law2_ScGeom_FrictPhys_Basic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* ncb){
	Real dt = Omega::instance().getTimeStep();
	int id1 = contact->getId1(), id2 = contact->getId2();

	ScGeom*    currentContactGeometry= static_cast<ScGeom*>(ig.get());
	FrictPhys* currentContactPhysics = static_cast<FrictPhys*>(ip.get());
	if(currentContactGeometry->penetrationDepth <0){
		if (neverErase) {
			currentContactPhysics->shearForce = Vector3r::Zero();
			currentContactPhysics->normalForce = Vector3r::Zero();}
		else 	ncb->interactions->requestErase(id1,id2);
		return;}
	State* de1 = Body::byId(id1,ncb)->state.get();
	State* de2 = Body::byId(id2,ncb)->state.get();
	Vector3r& shearForce = currentContactPhysics->shearForce;
	Real un=currentContactGeometry->penetrationDepth;
	TRVAR3(currentContactGeometry->penetrationDepth,de1->se3.position,de2->se3.position);
	currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal;
	
	if (!traceEnergy){//Update force but don't compute energy terms
		if(useShear){
			currentContactGeometry->updateShear(de1,de2,dt);
			shearForce=currentContactPhysics->ks*currentContactGeometry->shear;
		} else {
			currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);}
		// PFC3d SlipModel, is using friction angle. CoulombCriterion
		Real maxFs = currentContactPhysics->normalForce.squaredNorm()*
			std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
		if( shearForce.squaredNorm() > maxFs ){
			Real ratio = Mathr::Sqrt(maxFs) / shearForce.norm();
			shearForce *= ratio;
			if(useShear) currentContactGeometry->shear*=ratio;}
	} else {//almost the same with 2 additional Vector3r instanciated for energy tracing, duplicated block to make sure there is no cost for the instanciation of the vectors when traceEnergy==false
		if(useShear) throw ("energy tracing not defined with useShear==true");
		/*{
			currentContactGeometry->updateShear(de1,de2,dt);
			shearForce=currentContactPhysics->ks*currentContactGeometry->shear;//FIXME : energy terms if useShear?
		} else {*/
		Vector3r shearDisp = currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);
 		//}
		Real maxFs = currentContactPhysics->normalForce.squaredNorm()*
			std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
		if( shearForce.squaredNorm() > maxFs ){
			Real ratio = Mathr::Sqrt(maxFs) / shearForce.norm();
			Vector3r trialForce=shearForce;//store prev force for definition of plastic slip
			//define the plastic work input and increment the total plastic energy dissipated
			shearForce *= ratio;
			plasticDissipation +=
			(/*shearDisp+*/(1/currentContactPhysics->ks)*(trialForce-shearForce))//plastic disp.
			.dot(shearForce);//active force			
			//if(useShear) currentContactGeometry->shear*=ratio;
		}
	}	
	applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position, ncb);
	currentContactPhysics->prevNormal = currentContactGeometry->normal;
}

// same as elasticContactLaw, but using Dem3DofGeom
void Law2_Dem3DofGeom_FrictPhys_Basic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene*){
	Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(ig.get());
	FrictPhys* phys=static_cast<FrictPhys*>(ip.get());
	Real displN=geom->displacementN();
	if(displN>0){ scene->interactions->requestErase(contact->getId1(),contact->getId2()); return; }
	phys->normalForce=phys->kn*displN*geom->normal;
	Real maxFsSq=phys->normalForce.squaredNorm()*pow(phys->tangensOfFrictionAngle,2);
	Vector3r trialFs=phys->ks*geom->displacementT();
	if(trialFs.squaredNorm()>maxFsSq){ geom->slipToDisplacementTMax(sqrt(maxFsSq)/phys->ks); trialFs*=sqrt(maxFsSq/(trialFs.squaredNorm()));}
	phys->shearForce=trialFs;
	applyForceAtContactPoint(phys->normalForce+trialFs,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,scene);
}

// same as elasticContactLaw, but using Dem3DofGeom
void Law2_Dem6DofGeom_FrictPhys_Beam::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* scene){
	// normal & shear forces
	Dem6DofGeom* geom=static_cast<Dem6DofGeom*>(ig.get());
	FrictPhys* phys=static_cast<FrictPhys*>(ip.get());
	Real displN=geom->displacementN();
	phys->normalForce=phys->kn*displN*geom->normal;
	phys->shearForce=phys->ks*geom->displacementT();
	applyForceAtContactPoint(phys->normalForce+phys->shearForce,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,scene);
	// bend&twist:
	Vector3r bend; Real twist;
	geom->bendTwistAbs(bend,twist);
	Vector3r tt=bend*phys->kn+geom->normal*twist*phys->kn;
	cerr<<twist<<";"<<bend<<endl;
	scene->forces.addTorque(contact->getId1(),tt);
	scene->forces.addTorque(contact->getId2(),-tt);
}

/*************************************************************************
*  Copyright (C) 2004 by Olivier Galizzi                                 *
*  olivier.galizzi@xxxxxxx                                               *
*                                                                        *
*  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/GlobalEngine.hpp>

// only to see whether SCG_SHEAR is defined, may be removed in the future
#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/pkg-common/LawFunctor.hpp>
#include<yade/lib-base/openmp-accu.hpp>
#include <set>
#include <boost/tuple/tuple.hpp>


class Law2_ScGeom_FrictPhys_Basic: public LawFunctor{
	public:
		static Real Real0;
		OpenMPAccumulator<Real,&Law2_ScGeom_FrictPhys_Basic::Real0> plasticDissipation;
		virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene*);
		Real elasticEnergy ();
		Real getPlasticDissipation();
		void initPlasticDissipation(Real initVal=0);
	YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_FrictPhys_Basic,LawFunctor,"Law for linear compression, without cohesion and Mohr-Coulomb plasticity surface.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_Basic`, which uses :yref:`Dem3DofGeom`.",
		((int,sdecGroupMask,1,"Bitmask for allowing collision between particles :yref:`Body::groupMask`"))
		((bool,neverErase,false,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
		((bool,useShear,false,"Use ScGeom::updateShear rather than ScGeom::updateShearForce for shear force computation."))
		((bool,traceEnergy,false,"If true, increment plastic dissipation each time there is a plastic slip at contact. See :yref:`Law2_ScGeom_FrictPhys_Basic::plasticDissipation` and :yref:`Law2_ScGeom_FrictPhys_Basic::initPlasticDissipation`"))
		((Real,plasticW,0,"Register plastic energy for correct initialisation of the OpenMP accumulator. This value is NOT UPDATED "))
		,
		initPlasticDissipation(plasticW);
  		,
		.def("elasticEnergy",&Law2_ScGeom_FrictPhys_Basic::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts")
		.def("plasticDissipation",&Law2_ScGeom_FrictPhys_Basic::getPlasticDissipation,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_ScGeom_FrictPhys_Basic::traceEnergy` is true. |yupdate|")	.def("initPlasticDissipation",&Law2_ScGeom_FrictPhys_Basic::initPlasticDissipation,"Initialize cummulated plastic dissipation to a value (0 by default).")
	);
	FUNCTOR2D(ScGeom,FrictPhys);
	DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(Law2_ScGeom_FrictPhys_Basic);

/* Constitutive law for linear compression, no tension, and linear plasticity surface.

This class serves also as tutorial and is documented in detail at

	https://yade-dem.org/index.php/ConstitutiveLawHowto
*/
class Law2_Dem3DofGeom_FrictPhys_Basic: public LawFunctor{
	public:
		virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene*);
		FUNCTOR2D(Dem3DofGeom,FrictPhys);
		YADE_CLASS_BASE_DOC(Law2_Dem3DofGeom_FrictPhys_Basic,LawFunctor,"Constitutive law for linear compression, no tension, and linear plasticity surface.\n\nThis class serves also as tutorial and is documented in detail at https://yade-dem.org/index.php/ConstitutiveLawHowto.";);
};
REGISTER_SERIALIZABLE(Law2_Dem3DofGeom_FrictPhys_Basic);

/* Class for demonstrating beam-like behavior of the contact (normal, shear, bend and twist) */
class Law2_Dem6DofGeom_FrictPhys_Beam: public LawFunctor{
	public:
		virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene*);
	FUNCTOR2D(Dem6DofGeom,FrictPhys);
	YADE_CLASS_BASE_DOC(Law2_Dem6DofGeom_FrictPhys_Beam,LawFunctor,"Class for demonstrating beam-like behavior of contact (normal, shear, bend and twist) [broken][experimental]");
};
REGISTER_SERIALIZABLE(Law2_Dem6DofGeom_FrictPhys_Beam);

class ElasticContactLaw : public GlobalEngine{
		shared_ptr<Law2_ScGeom_FrictPhys_Basic> functor;
	public :
		void action();
	YADE_CLASS_BASE_DOC_ATTRS(ElasticContactLaw,GlobalEngine,"[DEPRECATED] Loop over interactions applying :yref:`Law2_ScGeom_FrictPhys_Basic` on all interactions.\n\n.. note::\n  Use :yref:`InteractionDispatchers` and :yref:`Law2_ScGeom_FrictPhys_Basic` instead of this class for performance reasons.",
		((int,sdecGroupMask,1,"Bitmask for allowing collision between particles :yref:`Body::groupMask`"))
		((bool,neverErase,false,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
		((bool,useShear,false,"Use :yref:`ScGeom`::updateShear rather than :yref:`ScGeom`::updateShearForce for shear force computation."))
	);
};

REGISTER_SERIALIZABLE(ElasticContactLaw);