yade-dev team mailing list archive

[gladky_anton@xxxxxxxxxxxxxxx]

Author: gladky_anton
Date: 2009-07-21 18:04:11 +0200 (Tue, 21 Jul 2009)
New Revision: 1880

Modified:
   trunk/pkg/dem/RockPM.cpp
   trunk/py/utils.py
   trunk/scripts/test/regular-sphere-pack.py
Log:
1. In utils.py alignedFacetBox has been changed to facetBox. It is now possible to create arbitrarily-aligned box composed of facets
2. Updated regular-sphere-pack.py
3. Some small changes on RockPM.cpp



Modified: trunk/pkg/dem/RockPM.cpp
===================================================================
--- trunk/pkg/dem/RockPM.cpp	2009-07-19 20:24:24 UTC (rev 1879)
+++ trunk/pkg/dem/RockPM.cpp	2009-07-21 16:04:11 UTC (rev 1880)
@@ -65,16 +65,19 @@
 		}
 	} else {
 		if (phys->isCohesive) {
-			phys->normalForce=phys->kn*displN*geom->normal;
 			if (displN>(phys->lengthMaxTension)) {
 				//LOG_WARN(displN<<"__TENSION!!!__");
 				phys->isCohesive = false;
 				rbp1->isDamaged=true;
 				rbp2->isDamaged=true;
+				return;
 			} else {
-				applyForceAtContactPoint(phys->normalForce,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,rootBody);
+				phys->normalForce=phys->kn*displN*geom->normal;
+				Real maxFsSq=phys->normalForce.SquaredLength()*pow(phys->tanFrictionAngle,2);
+				Vector3r trialFs=phys->ks*geom->displacementT();
+				if(trialFs.SquaredLength()>maxFsSq){ geom->slipToDisplacementTMax(sqrt(maxFsSq)); trialFs*=sqrt(maxFsSq/(trialFs.SquaredLength()));} 
+				applyForceAtContactPoint(phys->normalForce+trialFs,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,rootBody);
 			}
-			return;
 		} else {
 			rootBody->interactions->requestErase(contact->getId1(),contact->getId2());
 			return;
@@ -90,7 +93,7 @@
 	if(interaction->interactionPhysics) return; 
 
 	Dem3DofGeom* contGeom=YADE_CAST<Dem3DofGeom*>(interaction->interactionGeometry.get());
-	
+	Omega& OO=Omega::instance();
 	assert(contGeom);
 	//LOG_WARN(Omega::instance().getCurrentIteration());
 	const shared_ptr<RpmMat>& rpm1=YADE_PTR_CAST<RpmMat>(pp1);
@@ -117,9 +120,9 @@
 	
 	initDistance = contGeom->displacementN();
 
-	if ((rpm1->exampleNumber==rpm2->exampleNumber)&&(initDistance<(contPhys->lengthMaxTension))&&(initCohesive)) {
+	if ((rpm1->exampleNumber==rpm2->exampleNumber)&&(initDistance<(contPhys->lengthMaxTension))&&(initCohesive)&&(OO.getCurrentIteration()==0)) {
 		contPhys->isCohesive=true;
-		//LOG_WARN("InitDistance="<<initDistance);
+		//LOG_WARN("InitDistance="<<initDistance<<"  "<<rpm1->exampleNumber<<"  "<<rpm2->exampleNumber<<" "<<OO.getCurrentIteration());
 		//LOG_WARN("lengthMaxCompression="<<contPhys->lengthMaxCompression);
 		//LOG_WARN("lengthMaxTension="<<contPhys->lengthMaxTension);
 	}

Modified: trunk/py/utils.py
===================================================================
--- trunk/py/utils.py	2009-07-19 20:24:24 UTC (rev 1879)
+++ trunk/py/utils.py	2009-07-21 16:04:11 UTC (rev 1880)
@@ -117,16 +117,45 @@
 	b.mold.postProcessAttributes()
 	return b
 
-def alignedFacetBox(center,extents,wallMask=63,**kw):
-	"""Create axis-aligned box composed of facets, with given center and extents. wallMask determines which walls will be created,
+def facetBox(center,extents,orientation=[1,0,0,0],wallMask=63,**kw):
+	"""Create arbitrarily-aligned box composed of facets, with given center, extents and orientation. wallMask determines which walls will be created,
 	in the order -x (1), +x (2), -y (4), +y (8), -z (16), +z (32). The numbers are ANDed; the default 63 means to create all walls.
 	Remaining **kw arguments are passed to utils.facet. The facets are oriented outwards from the box."""
+	mn,mx=[-extents[i] for i in 0,1,2],[extents[i] for i in 0,1,2]
+	def doWall(a,b,c,d):
+		return [facet((a,b,c),**kw),facet((a,c,d),**kw)]
+	ret=[]
+	qTemp = Quaternion(Vector3(orientation[0],orientation[1],orientation[2]),orientation[3])
+	A=qTemp.Rotate(Vector3(mn[0],mn[1],mn[2]))+center
+	B=qTemp.Rotate(Vector3(mx[0],mn[1],mn[2]))+center
+	C=qTemp.Rotate(Vector3(mx[0],mx[1],mn[2]))+center
+	D=qTemp.Rotate(Vector3(mn[0],mx[1],mn[2]))+center
+	E=qTemp.Rotate(Vector3(mn[0],mn[1],mx[2]))+center
+	F=qTemp.Rotate(Vector3(mx[0],mn[1],mx[2]))+center
+	G=qTemp.Rotate(Vector3(mx[0],mx[1],mx[2]))+center
+	H=qTemp.Rotate(Vector3(mn[0],mx[1],mx[2]))+center
+	if wallMask&1:  ret+=doWall(A,D,H,E)
+	if wallMask&2:  ret+=doWall(B,C,G,F)
+	if wallMask&4:  ret+=doWall(A,B,F,E)
+	if wallMask&8:  ret+=doWall(D,H,G,C)
+	if wallMask&16: ret+=doWall(A,D,C,B)
+	if wallMask&32: ret+=doWall(E,F,G,H)
+	return ret
+	
+	'''
 	mn,mx=[center[i]-extents[i] for i in 0,1,2],[center[i]+extents[i] for i in 0,1,2]
 	def doWall(a,b,c,d):
 		return [facet((a,b,c),**kw),facet((a,c,d),**kw)]
 	ret=[]
-	A,B,C,D=(mn[0],mn[1],mn[2]),(mx[0],mn[1],mn[2]),(mx[0],mx[1],mn[2]),(mn[0],mx[1],mn[2])
-	E,F,G,H=(mn[0],mn[1],mx[2]),(mx[0],mn[1],mx[2]),(mx[0],mx[1],mx[2]),(mn[0],mx[1],mx[2])
+	qTemp = Quaternion(Vector3(orientation[0],orientation[1],orientation[2]),orientation[3])
+	A=qTemp.Rotate(Vector3(mn[0],mn[1],mn[2]))
+	B=qTemp.Rotate(Vector3(mx[0],mn[1],mn[2]))
+	C=qTemp.Rotate(Vector3(mx[0],mx[1],mn[2]))
+	D=qTemp.Rotate(Vector3(mn[0],mx[1],mn[2]))
+	E=qTemp.Rotate(Vector3(mn[0],mn[1],mx[2]))
+	F=qTemp.Rotate(Vector3(mx[0],mn[1],mx[2]))
+	G=qTemp.Rotate(Vector3(mx[0],mx[1],mx[2]))
+	H=qTemp.Rotate(Vector3(mn[0],mx[1],mx[2]))
 	if wallMask&1:  ret+=doWall(A,D,H,E)
 	if wallMask&2:  ret+=doWall(B,C,G,F)
 	if wallMask&4:  ret+=doWall(A,B,F,E)
@@ -134,8 +163,9 @@
 	if wallMask&16: ret+=doWall(A,D,C,B)
 	if wallMask&32: ret+=doWall(E,F,G,H)
 	return ret
+	'''
+	
 
-
 def aabbWalls(extrema=None,thickness=None,oversizeFactor=1.5,**kw):
 	"""return 6 walls that will wrap existing packing;
 	extrema are extremal points of the AABB of the packing (will be calculated if not specified)

Modified: trunk/scripts/test/regular-sphere-pack.py
===================================================================
--- trunk/scripts/test/regular-sphere-pack.py	2009-07-19 20:24:24 UTC (rev 1879)
+++ trunk/scripts/test/regular-sphere-pack.py	2009-07-21 16:04:11 UTC (rev 1880)
@@ -38,6 +38,19 @@
 	pack.regularOrtho(pack.inCylinder((-2,0,2),(-5,0,4),1),radius=rad,gap=gap,color=(0,0,1),**kw) # left hand
 	]: O.bodies.appendClumped(part)
 
+
+kwBoxes={'frictionAngle':0.5,'color':[1,0,0],'wire':False,'dynamic':False}
+q1 = Quaternion(Vector3(0,0,1),(3.14159/3))
+o1,o_angl = q1.ToAxisAngle()
+O.bodies.append(utils.facetBox((12,0,-6+0.9),(1,0.7,0.9),(o1[0],o1[1],o1[2],o_angl),**kwBoxes))
+q1 = Quaternion(Vector3(0,0,1),(3.14159/2))
+o1,o_angl = q1.ToAxisAngle()
+O.bodies.append(utils.facetBox((0,12,-6+0.9),(1,0.7,0.9),(o1[0],o1[1],o1[2],o_angl),**kwBoxes))
+q1 = Quaternion(Vector3(0,0,1),(3.14159))
+o1,o_angl = q1.ToAxisAngle()
+O.bodies.append(utils.facetBox((-12,-12,-6+0.9),(1,0.7,0.9),(o1[0],o1[1],o1[2],o_angl),**kwBoxes))
+
+
 try:
 	from yade import qt
 	qt.Controller()