yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3847
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Tue 2016-04-19 03:32:36 +0200
message:
  Modified MortarMat
  
  code restructuralization
  deleted some variables
  added example
added:
  examples/mortar/matModel2.py
modified:
  pkg/dem/MortarMat.cpp
  pkg/dem/MortarMat.hpp


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== added file 'examples/mortar/matModel2.py'
--- examples/mortar/matModel2.py	1970-01-01 00:00:00 +0000
+++ examples/mortar/matModel2.py	2016-04-19 01:32:36 +0000
@@ -0,0 +1,54 @@
+from yade import plot
+
+dt = 1
+young=1e9
+
+x,y,z = .2,.5,.7
+
+mortar = MortarMat(young=young)
+polyMat = PolyhedraMat(young=young/z)
+for mat in (mortar,polyMat):
+	O.materials.append(mat)
+
+O.dt = dt
+bs = b1,b2 = [polyhedron(((-x,-y,-z),(+x,-y,-z),(-x,+y,-z),(+x,+y,-z),(-x,-y,+z),(+x,-y,+z),(-x,+y,+z),(+x,+y,+z)),material=mortar) for i in (0,1)]
+b2.state.pos = b2.state.refPos = (0,0,2*z)
+for b in bs:
+	b.state.blockedDOFs = 'xyzXYZ'
+O.bodies.append(bs)
+
+def plotAddData():
+	i = O.interactions[0,1]
+	if i.phys:
+		plot.addData(
+			fn = i.phys.normalForce.norm(),
+			dspl = O.bodies[1].state.displ().norm(),
+		)
+
+#
+factor=1.1
+O.engines=[
+	ForceResetter(),
+	InsertionSortCollider([Bo1_Polyhedra_Aabb(aabbEnlargeFactor=factor,label='bo1')]),
+	InteractionLoop(
+		[Ig2_Polyhedra_Polyhedra_PolyhedraGeomOrScGeom(label='ig2')], 
+		[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(),Ip2_MortarMat_MortarMat_MortarPhys()],
+		[Law2_PolyhedraGeom_PolyhedraPhys_Volumetric(),Law2_ScGeom_MortarPhys_Lourenco()],
+	),
+	NewtonIntegrator(),
+	PyRunner(iterPeriod=1,command='plotAddData()'),
+	PyRunner(iterPeriod=1,command='print O.interactions[0,1].phys'),
+]
+ig2.ig2scGeom.interactionDetectionFactor = factor
+O.step()
+phys = O.interactions[0,1].phys
+phys.kn = young*2*x*2*y/(2*z)
+ig2.ig2scGeom.interactionDetectionFactor = bo1.aabbEnlargeFactor = 1
+ig2.createScGeom=False
+b2.state.vel = (0,0,-z*.002)
+for b in O.bodies:
+	b.mat = polyMat
+
+plot.plots = {'dspl':'fn'}
+plot.plot()
+O.run(30)

=== modified file 'pkg/dem/MortarMat.cpp'
--- pkg/dem/MortarMat.cpp	2016-04-13 16:45:33 +0000
+++ pkg/dem/MortarMat.cpp	2016-04-19 01:32:36 +0000
@@ -8,11 +8,14 @@
 CREATE_LOGGER(Ip2_MortarMat_MortarMat_MortarPhys);
 void Ip2_MortarMat_MortarMat_MortarPhys::go(const shared_ptr<Material>& material1, const shared_ptr<Material>& material2, const shared_ptr<Interaction>& interaction){
 	if (interaction->phys) return;
-	if (scene->iter >= cohesiveThresholdIter) return;
+	if (scene->iter >= cohesiveThresholdIter) {
+		LOG_ERROR("MortarMat not implemented for non-cohesive contacts");
+	}
 	shared_ptr<MortarPhys> phys(new MortarPhys());
 	interaction->phys = phys;
 	MortarMat* mat1 = YADE_CAST<MortarMat*>(material1.get());
 	MortarMat* mat2 = YADE_CAST<MortarMat*>(material2.get());
+	GenericSpheresContact* geom = YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
 
 	if (mat1->id>=0 && mat1->id == mat2->id) {
 		#define _CPATTR(a) phys->a=mat1->a
@@ -20,6 +23,7 @@
 			_CPATTR(compressiveStrength);
 			_CPATTR(cohesion);
 			_CPATTR(ellAspect);
+			_CPATTR(neverDamage);
 		#undef _CPATTR
 		phys->tangensOfFrictionAngle = std::tan(mat1->frictionAngle);
 	} else {
@@ -32,17 +36,36 @@
 			_AVGATTR(ellAspect);
 		#undef _AVGATTR
 		#undef _MINATTR
+		phys->neverDamage = mat1->neverDamage || mat2->neverDamage;
 		phys->tangensOfFrictionAngle = std::tan(.5*(mat1->frictionAngle+mat2->frictionAngle));
-		// E, G, kn, ks, crosssection, refPD, refLength to be computed in Law2
 	}
+	//
+	const Real& r1 = geom->refR1;
+	const Real& r2 = geom->refR2;
+	Real minRad = r1 <= 0 ? r2 : r2 <= 0 ? r1 : std::min(r1,r2);
+	phys->crossSection = std::pow(minRad,2);
+	const Real& E1 = mat1->young;
+	const Real& E2 = mat2->young;
+	const Real& n1 = mat1->poisson;
+	const Real& n2 = mat2->poisson;
+	phys->kn = 2*E1*r1*E2*r2/(E1*r1+E2*r2);
+	phys->ks = 2*E1*r1*n1*E2*r2*n2/(E1*r1*n1+E2*r2*n2);
 }
 
 
 
 
 CREATE_LOGGER(MortarPhys);
+
 MortarPhys::~MortarPhys(){};
 
+bool MortarPhys::failureCondition(Real sigmaN, Real sigmaT) {
+	bool cond1 = sigmaN - tensileStrength > 0;
+	bool cond2 = sigmaT + sigmaN*tangensOfFrictionAngle - cohesion > 0;
+	bool cond3 = std::pow(sigmaN,2) + std::pow(ellAspect*sigmaT,2) - std::pow(compressiveStrength,2) > 0;
+	return cond1 || cond2 || cond3;
+}
+
 
 
 
@@ -60,68 +83,38 @@
 	const shared_ptr<Body> b1 = Body::byId(id1,scene);
 	const shared_ptr<Body> b2 = Body::byId(id2,scene);
 
-	/* just the first time */
-	if (interaction->isFresh(scene)) {
-		const Vector3r& pos1 = b1->state->pos;
-		const Vector3r& pos2 = b2->state->pos;
-		const Real& r1 = geom->refR1;
-		const Real& r2 = geom->refR2;
-		Vector3r shift2 = scene->isPeriodic? Vector3r(scene->cell->hSize*interaction->cellDist.cast<Real>()) : Vector3r::Zero();
-		phys->refLength = (pos2 - pos1 + shift2).norm();
-		Real minRad = r1 <= 0 ? r2 : r2 <= 0 ? r1 : std::min(r1,r2);
-		phys->crossSection = std::pow(minRad,2);
-		phys->refPD = geom->refR1 + geom->refR2 - phys->refLength;
-		const shared_ptr<MortarMat> mat1 = YADE_PTR_CAST<MortarMat>(b1->material);
-		const shared_ptr<MortarMat> mat2 = YADE_PTR_CAST<MortarMat>(b2->material);
-		const Real& E1(mat1->young);
-		const Real& E2(mat2->young);
-		const Real& n1(mat1->poisson);
-		const Real& n2(mat2->poisson);
-		phys->kn = 2*E1*r1*E2*r2/(E1*r1+E2*r2);
-		phys->ks = 2*E1*r1*n1*E2*r2*n2/(E1*r1*n1+E2*r2*n2);
-		phys->E = phys->kn * phys->refLength / phys->crossSection;
-		phys->G = phys->ks * phys->refLength / phys->crossSection;
-	}
-	
 	/* shorthands */
-	Real& epsN(phys->epsN);
-	Vector3r& epsT(phys->epsT);
-	const Real& E(phys->E); \
-	const Real& G(phys->G);
 	const Real& crossSection(phys->crossSection);
 	Real& sigmaN(phys->sigmaN);
 	Vector3r& sigmaT(phys->sigmaT);
+	Real& kn(phys->kn);
+	Real& ks(phys->ks);
+	Vector3r& normalForce(phys->normalForce);
+	Vector3r& shearForce(phys->shearForce);
 
-	epsN = - (-phys->refPD + geom->penetrationDepth) / phys->refLength;
-	geom->rotate(epsT);
-	epsT += geom->shearIncrement() / phys->refLength;
 
 	/* constitutive law */
-	sigmaN = E*epsN;
-	sigmaT = G*epsT;
+	normalForce = kn*geom->penetrationDepth*geom->normal;
+	geom->rotate(shearForce);
+	shearForce -= ks*geom->shearIncrement();
+	sigmaN = - normalForce.dot(geom->normal) / crossSection;
+	sigmaT = - shearForce / crossSection;
 
-	Real st = sigmaT.norm();
-	bool cond1 = sigmaN - phys->tensileStrength > 0;
-	bool cond2 = st + sigmaN*phys->tangensOfFrictionAngle - phys->cohesion > 0;
-	bool cond3 = std::pow(sigmaN,2) + std::pow(phys->ellAspect*st,2) - std::pow(phys->compressiveStrength,2) > 0;
-	if (cond1 || cond2 || cond3) {
+	if (!phys->neverDamage && phys->failureCondition(sigmaN,sigmaT.norm())) {
 		return false;
 	}
    
-	phys->normalForce = -sigmaN*crossSection*geom->normal;
-	phys->shearForce = -sigmaT*crossSection;
-
 	State* s1 = b1->state.get();
 	State* s2 = b2->state.get();	
 
-	Vector3r f = -phys->normalForce - phys->shearForce;
+	Vector3r f = - normalForce - shearForce;
 	if (!scene->isPeriodic) {
 		applyForceAtContactPoint(f, geom->contactPoint , id1, s1->se3.position, id2, s2->se3.position);
 	} else {
 		scene->forces.addForce(id1,f);
 		scene->forces.addForce(id2,-f);
-		scene->forces.addTorque(id1,(geom->radius1+.5*(phys->refPD-geom->penetrationDepth))*geom->normal.cross(f));
-		scene->forces.addTorque(id2,(geom->radius2+.5*(phys->refPD-geom->penetrationDepth))*geom->normal.cross(f));
+		scene->forces.addTorque(id1,(geom->radius1-.5*(geom->penetrationDepth))*geom->normal.cross(f));
+		scene->forces.addTorque(id2,(geom->radius2-.5*(geom->penetrationDepth))*geom->normal.cross(f));
 	}
 	return true;
 }

=== modified file 'pkg/dem/MortarMat.hpp'
--- pkg/dem/MortarMat.hpp	2016-04-13 16:45:33 +0000
+++ pkg/dem/MortarMat.hpp	2016-04-19 01:32:36 +0000
@@ -26,6 +26,7 @@
 			((Real,compressiveStrength,10e6,,"compressiveStrength [Pa]"))
 			((Real,cohesion,1e6,,"cohesion [Pa]"))
 			((Real,ellAspect,3,,"aspect ratio of elliptical 'cap'. Value >1 means the ellipse is longer along normal stress axis."))
+			((bool,neverDamage,false,,"If true, interactions remain elastic regardless stresses"))
 			,
 			createIndex();
 		);
@@ -40,27 +41,23 @@
 
 class MortarPhys: public FrictPhys {
 	public:
-		Real epsN, sigmaN;
-		Vector3r epsT, sigmaT;
+		Real sigmaN;
+		Vector3r sigmaT;
 		virtual ~MortarPhys();
+		bool failureCondition(Real sigmaN, Real sigmaT);
 		YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(MortarPhys,FrictPhys,"IPhys class containing parameters of MortarMat. Used by Law2_ScGeom_MortarPhys_Lourenco.",
 			((Real,tensileStrength,NaN,,"tensileStrength [Pa]"))
 			((Real,compressiveStrength,NaN,,"compressiveStrength [Pa]"))
 			((Real,cohesion,NaN,,"cohesion [Pa]"))
 			((Real,ellAspect,NaN,,"aspect ratio of elliptical 'cap'. Value >1 means the ellipse is longer along normal stress axis."))
 			((Real,crossSection,NaN,,"Crosssection of interaction"))
-			((Real,refLength,NaN,,"Initial length of interaction [m]"))
-			((Real,refPD,NaN,,"Reference penetration depth [m]"))
-			((Real,E,NaN,,"interaction Young's modulus [Pa]"))
-			((Real,G,NaN,,"interaction shear modulus [Pa]"))
+			((bool,neverDamage,false,,"If true, interactions remain elastic regardless stresses"))
 			, // ctors
-			epsT = Vector3r::Zero();
 			createIndex();
 			,
-			.def_readonly("epsN",&MortarPhys::epsN,"Current normal strain |yupdate|")
-			.def_readonly("epsT",&MortarPhys::epsT,"Current shear strain |yupdate|")
 			.def_readonly("sigmaN",&MortarPhys::sigmaN,"Current normal stress |yupdate|")
 			.def_readonly("sigmaT",&MortarPhys::sigmaT,"Current shear stress |yupdate|")
+			.def("failureCondition",&MortarPhys::failureCondition,"Failure condition from normal stress and norm of shear stress (false=elastic, true=damaged)")
 		);
 		DECLARE_LOGGER;
 		REGISTER_CLASS_INDEX(MortarPhys,NormShearPhys);