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()
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
