yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10196
[Branch ~yade-pkg/yade/git-trunk] Rev 3742: Quite important changes in JCFpm code
------------------------------------------------------------
revno: 3742
committer: Jerome Duriez <jerome.duriez@xxxxxxxxxxxxxxx>
timestamp: Mon 2013-11-04 15:20:24 +0100
message:
Quite important changes in JCFpm code
- formal changes : all mechanical parameters for joint interactions, previously defined as parameters of Ip2..., are now stored in JCFpmMat. Sounds more logical for me and allows setting different joint parameters to different bodies.
- not formal change : the local stifnesses for joint interactions are now computed according to some surface relative to the in-contact spheres (instead of their radii), moreover the corresponding mechanical parameter. This parameter has the unit of Pa/m, rather than Pa, which may sound more logical for interface media.
modified:
examples/jointedCohesiveFrictionalPM/gravityBis.py
examples/jointedCohesiveFrictionalPM/identifBis.py
examples/jointedCohesiveFrictionalPM/identificationSpheresOnJoint.py
pkg/dem/JointedCohesiveFrictionalPM.cpp
pkg/dem/JointedCohesiveFrictionalPM.hpp
py/export.py
--
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
=== modified file 'examples/jointedCohesiveFrictionalPM/gravityBis.py'
--- examples/jointedCohesiveFrictionalPM/gravityBis.py 2013-10-14 08:28:54 +0000
+++ examples/jointedCohesiveFrictionalPM/gravityBis.py 2013-11-04 14:20:24 +0000
@@ -13,7 +13,7 @@
# the packing of spheres :
-def mat(): return JCFpmMat(type=1,young=1e8,frictionAngle=radians(30),density=3000)
+def mat(): return JCFpmMat(type=1,young=1e8,poisson=0.3,frictionAngle=radians(30),density=3000)
nSpheres = 3000.0
poros=0.13
rMeanSpheres = dimModele * pow(3.0/4.0*(1-poros)/(pi*nSpheres),1.0/3.0)
=== modified file 'examples/jointedCohesiveFrictionalPM/identifBis.py'
--- examples/jointedCohesiveFrictionalPM/identifBis.py 2013-10-04 09:47:21 +0000
+++ examples/jointedCohesiveFrictionalPM/identifBis.py 2013-11-04 14:20:24 +0000
@@ -13,7 +13,7 @@
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='is2aabb'),Bo1_Facet_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='ss2d3dg'),Ig2_Facet_Sphere_ScGeom()],
- [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,alpha=1,tensileStrength=1e6,cohesion=1e6,jointNormalStiffness=1,jointShearStiffness=1,jointTensileStrength=1e6,jointCohesion=1e6,jointFrictionAngle=1,label='interactionPhys')],
+ [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True,label='interactionLaw')]
),
NewtonIntegrator(damping=1)
@@ -47,31 +47,31 @@
nRef=normalRef/(normalRef.norm()) ## normalizes normalRef
normalFacetSphere=i.geom.normal # geom.normal is oriented from id1 to id2 -> normalFacetSphere from facet (i.id1) to sphere (i.id2)
- if O.bodies[i.id2].mat.onJoint==False : ## particles has not yet been identified as belonging to a joint plane
- O.bodies[i.id2].mat.onJoint=True
- O.bodies[i.id2].mat.joint=1
+ if O.bodies[i.id2].state.onJoint==False : ## particles has not yet been identified as belonging to a joint plane
+ O.bodies[i.id2].state.onJoint=True
+ O.bodies[i.id2].state.joint=1
O.bodies[i.id2].shape.color=jointcolor1
if nRef.dot(normalFacetSphere)>=0 :
- O.bodies[i.id2].mat.jointNormal1=nRef
+ O.bodies[i.id2].state.jointNormal1=nRef
elif nRef.dot(normalFacetSphere)<0 :
- O.bodies[i.id2].mat.jointNormal1=-nRef
- elif O.bodies[i.id2].mat.onJoint==True : ## particles has already been identified as belonging to, at least, 1 facet
- if O.bodies[i.id2].mat.joint==1 and ((O.bodies[i.id2].mat.jointNormal1.cross(nRef)).norm()>0.05) : ## particles has already been identified as belonging to only 1 facet
- O.bodies[i.id2].mat.joint=2
+ O.bodies[i.id2].state.jointNormal1=-nRef
+ elif O.bodies[i.id2].state.onJoint==True : ## particles has already been identified as belonging to, at least, 1 facet
+ if O.bodies[i.id2].state.joint==1 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) : ## particles has already been identified as belonging to only 1 facet
+ O.bodies[i.id2].state.joint=2
O.bodies[i.id2].shape.color=jointcolor2
if nRef.dot(normalFacetSphere)>=0 :
- O.bodies[i.id2].mat.jointNormal2=nRef
+ O.bodies[i.id2].state.jointNormal2=nRef
elif nRef.dot(normalFacetSphere)<0 :
- O.bodies[i.id2].mat.jointNormal2=-nRef
- elif O.bodies[i.id2].mat.joint==2 and ((O.bodies[i.id2].mat.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].mat.jointNormal2.cross(nRef)).norm()>0.05): ## particles has already been identified as belonging to more than 1 facet
- O.bodies[i.id2].mat.joint=3
+ O.bodies[i.id2].state.jointNormal2=-nRef
+ elif O.bodies[i.id2].state.joint==2 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].state.jointNormal2.cross(nRef)).norm()>0.05): ## particles has already been identified as belonging to more than 1 facet
+ O.bodies[i.id2].state.joint=3
O.bodies[i.id2].shape.color=jointcolor3
if nRef.dot(normalFacetSphere)>=0 :
- O.bodies[i.id2].mat.jointNormal3=nRef
+ O.bodies[i.id2].state.jointNormal3=nRef
elif nRef.dot(normalFacetSphere)<0 :
- O.bodies[i.id2].mat.jointNormal3=-nRef
- elif O.bodies[i.id2].mat.joint==3 and ((O.bodies[i.id2].mat.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].mat.jointNormal2.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].mat.jointNormal3.cross(nRef)).norm()>0.05):
- O.bodies[i.id2].mat.joint=4
+ O.bodies[i.id2].state.jointNormal3=-nRef
+ elif O.bodies[i.id2].state.joint==3 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].state.jointNormal2.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].state.jointNormal3.cross(nRef)).norm()>0.05):
+ O.bodies[i.id2].state.joint=4
O.bodies[i.id2].shape.color=jointcolor5
#### second step -> find spheres interacting with spheres on facet (could be executed in the same timestep as step 1?)
@@ -82,12 +82,12 @@
vertices=O.bodies[j.id1].shape.vertices
normalRef=vertices[0].cross(vertices[1]) # defines the normal to the facet normalRef
nRef=normalRef/(normalRef.norm()) ## normalizes normalRef
- if ((O.bodies[j.id2].mat.jointNormal1.cross(nRef)).norm()<0.05) :
- jointNormalRef=O.bodies[j.id2].mat.jointNormal1
- elif ((O.bodies[j.id2].mat.jointNormal2.cross(nRef)).norm()<0.05) :
- jointNormalRef=O.bodies[j.id2].mat.jointNormal2
- elif ((O.bodies[j.id2].mat.jointNormal3.cross(nRef)).norm()<0.05) :
- jointNormalRef=O.bodies[j.id2].mat.jointNormal3
+ if ((O.bodies[j.id2].state.jointNormal1.cross(nRef)).norm()<0.05) :
+ jointNormalRef=O.bodies[j.id2].state.jointNormal1
+ elif ((O.bodies[j.id2].state.jointNormal2.cross(nRef)).norm()<0.05) :
+ jointNormalRef=O.bodies[j.id2].state.jointNormal2
+ elif ((O.bodies[j.id2].state.jointNormal3.cross(nRef)).norm()<0.05) :
+ jointNormalRef=O.bodies[j.id2].state.jointNormal3
else : continue
facetCenter=O.bodies[j.id1].state.pos
#### seek for each sphere interacting with the identified sphere j.id2
@@ -100,51 +100,51 @@
sphOnF=n.id2
othSph=n.id1
facetSphereDir=(O.bodies[othSph].state.pos-facetCenter)
- if O.bodies[othSph].mat.onJoint==True :
- if O.bodies[othSph].mat.joint==3 and ((O.bodies[othSph].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].mat.jointNormal2.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].mat.jointNormal3.cross(jointNormalRef)).norm()>0.05):
- O.bodies[othSph].mat.joint=4
+ if O.bodies[othSph].state.onJoint==True :
+ if O.bodies[othSph].state.joint==3 and ((O.bodies[othSph].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].state.jointNormal2.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].state.jointNormal3.cross(jointNormalRef)).norm()>0.05):
+ O.bodies[othSph].state.joint=4
O.bodies[othSph].shape.color=jointcolor5
- elif O.bodies[othSph].mat.joint==2 and ((O.bodies[othSph].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].mat.jointNormal2.cross(jointNormalRef)).norm()>0.05):
- O.bodies[othSph].mat.joint=3
- if facetSphereDir.dot(jointNormalRef)>=0:
- O.bodies[othSph].mat.jointNormal3=jointNormalRef
- elif facetSphereDir.dot(jointNormalRef)<0:
- O.bodies[othSph].mat.jointNormal3=-jointNormalRef
- elif O.bodies[othSph].mat.joint==1 and ((O.bodies[othSph].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) :
- O.bodies[othSph].mat.joint=2
- if facetSphereDir.dot(jointNormalRef)>=0:
- O.bodies[othSph].mat.jointNormal2=jointNormalRef
- elif facetSphereDir.dot(jointNormalRef)<0:
- O.bodies[othSph].mat.jointNormal2=-jointNormalRef
- elif O.bodies[othSph].mat.onJoint==False :
- O.bodies[othSph].mat.onJoint=True
- O.bodies[othSph].mat.joint=1
+ elif O.bodies[othSph].state.joint==2 and ((O.bodies[othSph].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].state.jointNormal2.cross(jointNormalRef)).norm()>0.05):
+ O.bodies[othSph].state.joint=3
+ if facetSphereDir.dot(jointNormalRef)>=0:
+ O.bodies[othSph].state.jointNormal3=jointNormalRef
+ elif facetSphereDir.dot(jointNormalRef)<0:
+ O.bodies[othSph].state.jointNormal3=-jointNormalRef
+ elif O.bodies[othSph].state.joint==1 and ((O.bodies[othSph].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) :
+ O.bodies[othSph].state.joint=2
+ if facetSphereDir.dot(jointNormalRef)>=0:
+ O.bodies[othSph].state.jointNormal2=jointNormalRef
+ elif facetSphereDir.dot(jointNormalRef)<0:
+ O.bodies[othSph].state.jointNormal2=-jointNormalRef
+ elif O.bodies[othSph].state.onJoint==False :
+ O.bodies[othSph].state.onJoint=True
+ O.bodies[othSph].state.joint=1
O.bodies[othSph].shape.color=jointcolor4
if facetSphereDir.dot(jointNormalRef)>=0:
- O.bodies[othSph].mat.jointNormal1=jointNormalRef
+ O.bodies[othSph].state.jointNormal1=jointNormalRef
elif facetSphereDir.dot(jointNormalRef)<0:
- O.bodies[othSph].mat.jointNormal1=-jointNormalRef
+ O.bodies[othSph].state.jointNormal1=-jointNormalRef
#### for visualization:
#bj=0
#vert=(0.,1.,0.)
#hor=(0.,1.,1.)
#for o in O.bodies:
- #if o.mat.onJoint==True : # or o.shape.name=='Facet':
+ #if o.state.onJoint==True : # or o.shape.name=='Facet':
##if o.shape.name=='Facet':
##o.shape.wire=True
##o.state.pos+=(0,50,0)
##bj+=1
- #if o.mat.jointNormal1.dot(hor)>0 :
+ #if o.state.jointNormal1.dot(hor)>0 :
##o.state.pos+=(0,50,0)
#o.shape.color=jointcolor1
- #elif o.mat.jointNormal1.dot(hor)<0 :
+ #elif o.state.jointNormal1.dot(hor)<0 :
##o.state.pos+=(0,55,0)
#o.shape.color=jointcolor2
#if o.mat.type>2 :
#bj+=1
#o.shape.color=jointcolor5
- ##print o.mat.jointNormal.dot(vert)
+ ##print o.state.jointNormal.dot(vert)
##### to delete facets (should be OK to start the simulation after that!
@@ -154,8 +154,7 @@
O.resetTime()
O.interactions.clear()
-print ''
-print 'IdentificationSpheresOnJoint executed ! Spheres onJoint (and so on...) detected, facets deleted'
-print ''
+print '\n IdentificationSpheresOnJoint executed ! Spheres onJoint (and so on...) detected, facets deleted\n\n'
+
=== modified file 'examples/jointedCohesiveFrictionalPM/identificationSpheresOnJoint.py'
--- examples/jointedCohesiveFrictionalPM/identificationSpheresOnJoint.py 2013-10-04 09:43:53 +0000
+++ examples/jointedCohesiveFrictionalPM/identificationSpheresOnJoint.py 2013-11-04 14:20:24 +0000
@@ -81,31 +81,31 @@
nRef=normalRef/(normalRef.norm()) ## normalizes normalRef
normalFacetSphere=i.geom.normal # geom.normal is oriented from id1 to id2 -> normalFacetSphere from facet (i.id1) to sphere (i.id2)
- if O.bodies[i.id2].mat.onJoint==False : ## particles has not yet been identified as belonging to a joint plane
- O.bodies[i.id2].mat.onJoint=True
- O.bodies[i.id2].mat.joint=1
+ if O.bodies[i.id2].state.onJoint==False : ## particles has not yet been identified as belonging to a joint plane
+ O.bodies[i.id2].state.onJoint=True
+ O.bodies[i.id2].state.joint=1
O.bodies[i.id2].shape.color=jointcolor1
if nRef.dot(normalFacetSphere)>=0 :
- O.bodies[i.id2].mat.jointNormal1=nRef
+ O.bodies[i.id2].state.jointNormal1=nRef
elif nRef.dot(normalFacetSphere)<0 :
- O.bodies[i.id2].mat.jointNormal1=-nRef
- elif O.bodies[i.id2].mat.onJoint==True : ## particles has already been identified as belonging to, at least, 1 facet
- if O.bodies[i.id2].mat.joint==1 and ((O.bodies[i.id2].mat.jointNormal1.cross(nRef)).norm()>0.05) : ## particles has already been identified as belonging to only 1 facet
- O.bodies[i.id2].mat.joint=2
+ O.bodies[i.id2].state.jointNormal1=-nRef
+ elif O.bodies[i.id2].state.onJoint==True : ## particles has already been identified as belonging to, at least, 1 facet
+ if O.bodies[i.id2].state.joint==1 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) : ## particles has already been identified as belonging to only 1 facet
+ O.bodies[i.id2].state.joint=2
O.bodies[i.id2].shape.color=jointcolor2
if nRef.dot(normalFacetSphere)>=0 :
- O.bodies[i.id2].mat.jointNormal2=nRef
+ O.bodies[i.id2].state.jointNormal2=nRef
elif nRef.dot(normalFacetSphere)<0 :
- O.bodies[i.id2].mat.jointNormal2=-nRef
- elif O.bodies[i.id2].mat.joint==2 and ((O.bodies[i.id2].mat.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].mat.jointNormal2.cross(nRef)).norm()>0.05): ## particles has already been identified as belonging to more than 1 facet
- O.bodies[i.id2].mat.joint=3
+ O.bodies[i.id2].state.jointNormal2=-nRef
+ elif O.bodies[i.id2].state.joint==2 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].state.jointNormal2.cross(nRef)).norm()>0.05): ## particles has already been identified as belonging to more than 1 facet
+ O.bodies[i.id2].state.joint=3
O.bodies[i.id2].shape.color=jointcolor3
if nRef.dot(normalFacetSphere)>=0 :
- O.bodies[i.id2].mat.jointNormal3=nRef
+ O.bodies[i.id2].state.jointNormal3=nRef
elif nRef.dot(normalFacetSphere)<0 :
- O.bodies[i.id2].mat.jointNormal3=-nRef
- elif O.bodies[i.id2].mat.joint==3 and ((O.bodies[i.id2].mat.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].mat.jointNormal2.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].mat.jointNormal3.cross(nRef)).norm()>0.05):
- O.bodies[i.id2].mat.joint=4
+ O.bodies[i.id2].state.jointNormal3=-nRef
+ elif O.bodies[i.id2].state.joint==3 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].state.jointNormal2.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].state.jointNormal3.cross(nRef)).norm()>0.05):
+ O.bodies[i.id2].state.joint=4
O.bodies[i.id2].shape.color=jointcolor5
#### second step -> find spheres interacting with spheres on facet (could be executed in the same timestep as step 1?)
@@ -115,69 +115,69 @@
vertices=O.bodies[j.id1].shape.vertices
normalRef=vertices[0].cross(vertices[1]) # defines the normal to the facet normalRef
nRef=normalRef/(normalRef.norm()) ## normalizes normalRef
- if ((O.bodies[j.id2].mat.jointNormal1.cross(nRef)).norm()<0.05) :
- jointNormalRef=O.bodies[j.id2].mat.jointNormal1
- elif ((O.bodies[j.id2].mat.jointNormal2.cross(nRef)).norm()<0.05) :
- jointNormalRef=O.bodies[j.id2].mat.jointNormal2
- elif ((O.bodies[j.id2].mat.jointNormal3.cross(nRef)).norm()<0.05) :
- jointNormalRef=O.bodies[j.id2].mat.jointNormal3
+ if ((O.bodies[j.id2].state.jointNormal1.cross(nRef)).norm()<0.05) :
+ jointNormalRef=O.bodies[j.id2].state.jointNormal1
+ elif ((O.bodies[j.id2].state.jointNormal2.cross(nRef)).norm()<0.05) :
+ jointNormalRef=O.bodies[j.id2].state.jointNormal2
+ elif ((O.bodies[j.id2].state.jointNormal3.cross(nRef)).norm()<0.05) :
+ jointNormalRef=O.bodies[j.id2].state.jointNormal3
else : continue
facetCenter=O.bodies[j.id1].state.pos
#### seek for each sphere interacting with the identified sphere j.id2
for n in O.interactions.withBody(j.id2) :
if n.id1==j.id2 and isinstance(O.bodies[n.id2].shape,Sphere):
facetSphereDir=(O.bodies[n.id2].state.pos-facetCenter)
- if O.bodies[n.id2].mat.onJoint==True :
- if O.bodies[n.id2].mat.joint==3 and ((O.bodies[n.id2].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[n.id2].mat.jointNormal2.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[n.id2].mat.jointNormal3.cross(jointNormalRef)).norm()>0.05):
- O.bodies[n.id2].mat.joint=4
+ if O.bodies[n.id2].state.onJoint==True :
+ if O.bodies[n.id2].state.joint==3 and ((O.bodies[n.id2].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[n.id2].state.jointNormal2.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[n.id2].state.jointNormal3.cross(jointNormalRef)).norm()>0.05):
+ O.bodies[n.id2].state.joint=4
O.bodies[n.id2].shape.color=jointcolor5
- elif O.bodies[n.id2].mat.joint==2 and ((O.bodies[n.id2].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[n.id2].mat.jointNormal2.cross(jointNormalRef)).norm()>0.05):
- O.bodies[n.id2].mat.joint=3
- if facetSphereDir.dot(jointNormalRef)>=0:
- O.bodies[n.id2].mat.jointNormal3=jointNormalRef
- elif facetSphereDir.dot(jointNormalRef)<0:
- O.bodies[n.id2].mat.jointNormal3=-jointNormalRef
- elif O.bodies[n.id2].mat.joint==1 and ((O.bodies[n.id2].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) :
- O.bodies[n.id2].mat.joint=2
- if facetSphereDir.dot(jointNormalRef)>=0:
- O.bodies[n.id2].mat.jointNormal2=jointNormalRef
- elif facetSphereDir.dot(jointNormalRef)<0:
- O.bodies[n.id2].mat.jointNormal2=-jointNormalRef
- elif O.bodies[n.id2].mat.onJoint==False :
- O.bodies[n.id2].mat.onJoint=True
- O.bodies[n.id2].mat.joint=1
+ elif O.bodies[n.id2].state.joint==2 and ((O.bodies[n.id2].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[n.id2].state.jointNormal2.cross(jointNormalRef)).norm()>0.05):
+ O.bodies[n.id2].state.joint=3
+ if facetSphereDir.dot(jointNormalRef)>=0:
+ O.bodies[n.id2].state.jointNormal3=jointNormalRef
+ elif facetSphereDir.dot(jointNormalRef)<0:
+ O.bodies[n.id2].state.jointNormal3=-jointNormalRef
+ elif O.bodies[n.id2].state.joint==1 and ((O.bodies[n.id2].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) :
+ O.bodies[n.id2].state.joint=2
+ if facetSphereDir.dot(jointNormalRef)>=0:
+ O.bodies[n.id2].state.jointNormal2=jointNormalRef
+ elif facetSphereDir.dot(jointNormalRef)<0:
+ O.bodies[n.id2].state.jointNormal2=-jointNormalRef
+ elif O.bodies[n.id2].state.onJoint==False :
+ O.bodies[n.id2].state.onJoint=True
+ O.bodies[n.id2].state.joint=1
O.bodies[n.id2].shape.color=jointcolor4
if facetSphereDir.dot(jointNormalRef)>=0:
- O.bodies[n.id2].mat.jointNormal1=jointNormalRef
+ O.bodies[n.id2].state.jointNormal1=jointNormalRef
elif facetSphereDir.dot(jointNormalRef)<0:
- O.bodies[n.id2].mat.jointNormal1=-jointNormalRef
+ O.bodies[n.id2].state.jointNormal1=-jointNormalRef
elif n.id2==j.id2 and isinstance(O.bodies[n.id1].shape,Sphere):
facetSphereDir=(O.bodies[n.id1].state.pos-facetCenter)
- if O.bodies[n.id1].mat.onJoint==True :
- if O.bodies[n.id1].mat.joint==3 and ((O.bodies[n.id1].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[n.id1].mat.jointNormal2.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[n.id1].mat.jointNormal3.cross(jointNormalRef)).norm()>0.05):
- O.bodies[n.id1].mat.joint=4
+ if O.bodies[n.id1].state.onJoint==True :
+ if O.bodies[n.id1].state.joint==3 and ((O.bodies[n.id1].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[n.id1].state.jointNormal2.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[n.id1].state.jointNormal3.cross(jointNormalRef)).norm()>0.05):
+ O.bodies[n.id1].state.joint=4
O.bodies[n.id1].shape.color=jointcolor5
- elif O.bodies[n.id1].mat.joint==2 and ((O.bodies[n.id1].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[n.id1].mat.jointNormal2.cross(jointNormalRef)).norm()>0.05):
- O.bodies[n.id1].mat.joint=3
- if facetSphereDir.dot(jointNormalRef)>=0:
- O.bodies[n.id1].mat.jointNormal3=jointNormalRef
- elif facetSphereDir.dot(jointNormalRef)<0:
- O.bodies[n.id1].mat.jointNormal3=-jointNormalRef
- elif O.bodies[n.id1].mat.joint==1 and ((O.bodies[n.id1].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) :
- O.bodies[n.id1].mat.joint=2
- if facetSphereDir.dot(jointNormalRef)>=0:
- O.bodies[n.id1].mat.jointNormal2=jointNormalRef
- elif facetSphereDir.dot(jointNormalRef)<0:
- O.bodies[n.id1].mat.jointNormal2=-jointNormalRef
- elif O.bodies[n.id1].mat.onJoint==False :
- O.bodies[n.id1].mat.onJoint=True
- O.bodies[n.id1].mat.joint=1
+ elif O.bodies[n.id1].state.joint==2 and ((O.bodies[n.id1].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[n.id1].state.jointNormal2.cross(jointNormalRef)).norm()>0.05):
+ O.bodies[n.id1].state.joint=3
+ if facetSphereDir.dot(jointNormalRef)>=0:
+ O.bodies[n.id1].state.jointNormal3=jointNormalRef
+ elif facetSphereDir.dot(jointNormalRef)<0:
+ O.bodies[n.id1].state.jointNormal3=-jointNormalRef
+ elif O.bodies[n.id1].state.joint==1 and ((O.bodies[n.id1].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) :
+ O.bodies[n.id1].state.joint=2
+ if facetSphereDir.dot(jointNormalRef)>=0:
+ O.bodies[n.id1].state.jointNormal2=jointNormalRef
+ elif facetSphereDir.dot(jointNormalRef)<0:
+ O.bodies[n.id1].state.jointNormal2=-jointNormalRef
+ elif O.bodies[n.id1].state.onJoint==False :
+ O.bodies[n.id1].state.onJoint=True
+ O.bodies[n.id1].state.joint=1
O.bodies[n.id1].shape.color=jointcolor4
if facetSphereDir.dot(jointNormalRef)>=0:
- O.bodies[n.id1].mat.jointNormal1=jointNormalRef
+ O.bodies[n.id1].state.jointNormal1=jointNormalRef
elif facetSphereDir.dot(jointNormalRef)<0:
- O.bodies[n.id1].mat.jointNormal1=-jointNormalRef
+ O.bodies[n.id1].state.jointNormal1=-jointNormalRef
#### for visualization:
@@ -185,21 +185,21 @@
#vert=(0.,1.,0.)
#hor=(0.,1.,1.)
#for o in O.bodies:
- #if o.mat.onJoint==True : # or o.shape.name=='Facet':
+ #if o.state.onJoint==True : # or o.shape.name=='Facet':
##if o.shape.name=='Facet':
##o.shape.wire=True
##o.state.pos+=(0,50,0)
##bj+=1
- #if o.mat.jointNormal1.dot(hor)>0 :
+ #if o.state.jointNormal1.dot(hor)>0 :
##o.state.pos+=(0,50,0)
#o.shape.color=jointcolor1
- #elif o.mat.jointNormal1.dot(hor)<0 :
+ #elif o.state.jointNormal1.dot(hor)<0 :
##o.state.pos+=(0,55,0)
#o.shape.color=jointcolor2
#if o.mat.type>2 :
#bj+=1
#o.shape.color=jointcolor5
- ##print o.mat.jointNormal.dot(vert)
+ ##print o.state.jointNormal.dot(vert)
#### Save text file with informations on each sphere
=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.cpp'
--- pkg/dem/JointedCohesiveFrictionalPM.cpp 2013-11-04 14:03:39 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.cpp 2013-11-04 14:20:24 +0000
@@ -176,14 +176,22 @@
const shared_ptr<JCFpmMat>& yade1 = YADE_PTR_CAST<JCFpmMat>(b1);
const shared_ptr<JCFpmMat>& yade2 = YADE_PTR_CAST<JCFpmMat>(b2);
+ JCFpmState* st1=dynamic_cast<JCFpmState*>(Body::byId(interaction->getId1(),scene)->state.get());
+ JCFpmState* st2=dynamic_cast<JCFpmState*>(Body::byId(interaction->getId2(),scene)->state.get());
shared_ptr<JCFpmPhys> contactPhysics(new JCFpmPhys());
/* From material properties */
Real E1 = yade1->young;
Real E2 = yade2->young;
+ Real v1 = yade1->poisson;
+ Real v2 = yade2->poisson;
Real f1 = yade1->frictionAngle;
Real f2 = yade2->frictionAngle;
+ Real SigT1 = yade1->tensileStrength;
+ Real SigT2 = yade2->tensileStrength;
+ Real Coh1 = yade1->cohesion;
+ Real Coh2 = yade2->cohesion;
/* From interaction geometry */
Real R1= geom->radius1;
@@ -194,56 +202,82 @@
// frictional properties
contactPhysics->kn = 2.*E1*R1*E2*R2/(E1*R1+E2*R2);
- contactPhysics->ks = alpha*contactPhysics->kn;
- contactPhysics->frictionAngle = std::min(f1,f2);
- contactPhysics->tanFrictionAngle = std::tan(contactPhysics->frictionAngle);
+ contactPhysics->ks = 2.*E1*R1*v1*E2*R2*v2/(E1*R1*v1+E2*R2*v2);//alpha*contactPhysics->kn;
+ contactPhysics->tanFrictionAngle = std::tan(std::min(f1,f2));
// cohesive properties
///to set if the contact is cohesive or not
- if ((scene->iter < cohesiveTresholdIteration) && (tensileStrength>0 || cohesion>0) && (yade1->type == yade2->type)){
+ if ((scene->iter < cohesiveTresholdIteration) && (std::min(SigT1,SigT2)>0 || std::min(Coh1,Coh2)>0) && (yade1->type == yade2->type)){
contactPhysics->isCohesive=true;
- JCFpmState* st1=dynamic_cast<JCFpmState*>(Body::byId(interaction->getId1(),scene)->state.get());
st1->noIniLinks++;
- JCFpmState* st2=dynamic_cast<JCFpmState*>(Body::byId(interaction->getId2(),scene)->state.get());
st2->noIniLinks++;
}
if ( contactPhysics->isCohesive ) {
- contactPhysics->FnMax = tensileStrength*contactPhysics->crossSection;
- contactPhysics->FsMax = cohesion*contactPhysics->crossSection;
+ contactPhysics->FnMax = std::min(SigT1,SigT2)*contactPhysics->crossSection;
+ contactPhysics->FsMax = std::min(Coh1,Coh2)*contactPhysics->crossSection;
}
/// +++ Jointed interactions ->NOTE: geom->normal is oriented from 1 to 2 / jointNormal from plane to sphere
- if ( yade1->onJoint && yade2->onJoint ) {
-
- if ( (((yade1->jointNormal1.cross(yade2->jointNormal1)).norm()<0.1) && (yade1->jointNormal1.dot(yade2->jointNormal1)<0)) || (((yade1->jointNormal1.cross(yade2->jointNormal2)).norm()<0.1) && (yade1->jointNormal1.dot(yade2->jointNormal2)<0)) || (((yade1->jointNormal1.cross(yade2->jointNormal3)).norm()<0.1) && (yade1->jointNormal1.dot(yade2->jointNormal3)<0)) ) { contactPhysics->isOnJoint = true; contactPhysics->jointNormal = yade1->jointNormal1; }
- else if ( (((yade1->jointNormal2.cross(yade2->jointNormal1)).norm()<0.1) && (yade1->jointNormal2.dot(yade2->jointNormal1)<0)) || (((yade1->jointNormal2.cross(yade2->jointNormal2)).norm()<0.1) && (yade1->jointNormal2.dot(yade2->jointNormal2)<0)) || (((yade1->jointNormal2.cross(yade2->jointNormal3)).norm()<0.1) && (yade1->jointNormal2.dot(yade2->jointNormal3)<0)) ) { contactPhysics->isOnJoint = true; contactPhysics->jointNormal = yade1->jointNormal2; }
- else if ( (((yade1->jointNormal3.cross(yade2->jointNormal1)).norm()<0.1) && (yade1->jointNormal3.dot(yade2->jointNormal1)<0)) || (((yade1->jointNormal3.cross(yade2->jointNormal2)).norm()<0.1) && (yade1->jointNormal3.dot(yade2->jointNormal2)<0)) || (((yade1->jointNormal3.cross(yade2->jointNormal3)).norm()<0.1) && (yade1->jointNormal3.dot(yade2->jointNormal3)<0)) ) { contactPhysics->isOnJoint = true; contactPhysics->jointNormal = yade1->jointNormal3; }
- else if ( (yade1->joint>3 || yade2->joint>3) && ( ( ((yade1->jointNormal1.cross(yade2->jointNormal1)).norm()>0.1) && ((yade1->jointNormal1.cross(yade2->jointNormal2)).norm()>0.1) && ((yade1->jointNormal1.cross(yade2->jointNormal3)).norm()>0.1) ) || ( ((yade1->jointNormal2.cross(yade2->jointNormal1)).norm()>0.1) && ((yade1->jointNormal2.cross(yade2->jointNormal2)).norm()>0.1) && ((yade1->jointNormal2.cross(yade2->jointNormal3)).norm()>0.1) ) || ( ((yade1->jointNormal3.cross(yade2->jointNormal1)).norm()>0.1) && ((yade1->jointNormal3.cross(yade2->jointNormal2)).norm()>0.1) && ((yade1->jointNormal3.cross(yade2->jointNormal3)).norm()>0.1) ) ) ) { contactPhysics->isOnJoint = true; contactPhysics->more = true; contactPhysics->jointNormal = geom->normal; }
+ if ( st1->onJoint && st2->onJoint )
+ {
+ if ( (((st1->jointNormal1.cross(st2->jointNormal1)).norm()<0.1) && (st1->jointNormal1.dot(st2->jointNormal1)<0)) || (((st1->jointNormal1.cross(st2->jointNormal2)).norm()<0.1) && (st1->jointNormal1.dot(st2->jointNormal2)<0)) || (((st1->jointNormal1.cross(st2->jointNormal3)).norm()<0.1) && (st1->jointNormal1.dot(st2->jointNormal3)<0)) )
+ {
+ contactPhysics->isOnJoint = true;
+ contactPhysics->jointNormal = st1->jointNormal1;
+ }
+ else if ( (((st1->jointNormal2.cross(st2->jointNormal1)).norm()<0.1) && (st1->jointNormal2.dot(st2->jointNormal1)<0)) || (((st1->jointNormal2.cross(st2->jointNormal2)).norm()<0.1) && (st1->jointNormal2.dot(st2->jointNormal2)<0)) || (((st1->jointNormal2.cross(st2->jointNormal3)).norm()<0.1) && (st1->jointNormal2.dot(st2->jointNormal3)<0)) )
+ {
+ contactPhysics->isOnJoint = true;
+ contactPhysics->jointNormal = st1->jointNormal2;
+ }
+ else if ( (((st1->jointNormal3.cross(st2->jointNormal1)).norm()<0.1) && (st1->jointNormal3.dot(st2->jointNormal1)<0)) || (((st1->jointNormal3.cross(st2->jointNormal2)).norm()<0.1) && (st1->jointNormal3.dot(st2->jointNormal2)<0)) || (((st1->jointNormal3.cross(st2->jointNormal3)).norm()<0.1) && (st1->jointNormal3.dot(st2->jointNormal3)<0)) )
+ {
+ contactPhysics->isOnJoint = true;
+ contactPhysics->jointNormal = st1->jointNormal3;
+ }
+ else if ( (st1->joint>3 || st2->joint>3) && ( ( ((st1->jointNormal1.cross(st2->jointNormal1)).norm()>0.1) && ((st1->jointNormal1.cross(st2->jointNormal2)).norm()>0.1) && ((st1->jointNormal1.cross(st2->jointNormal3)).norm()>0.1) ) || ( ((st1->jointNormal2.cross(st2->jointNormal1)).norm()>0.1) && ((st1->jointNormal2.cross(st2->jointNormal2)).norm()>0.1) && ((st1->jointNormal2.cross(st2->jointNormal3)).norm()>0.1) ) || ( ((st1->jointNormal3.cross(st2->jointNormal1)).norm()>0.1) && ((st1->jointNormal3.cross(st2->jointNormal2)).norm()>0.1) && ((st1->jointNormal3.cross(st2->jointNormal3)).norm()>0.1) ) ) ) { contactPhysics->isOnJoint = true; contactPhysics->more = true; contactPhysics->jointNormal = geom->normal; }
}
///to specify joint properties
if ( contactPhysics->isOnJoint ) {
- contactPhysics->frictionAngle = (jointFrictionAngle>=0 ? jointFrictionAngle : std::min(f1,f2));
- contactPhysics->tanFrictionAngle = std::tan(contactPhysics->frictionAngle);
- contactPhysics->kn = jointNormalStiffness*2.*R1*R2/(R1+R2);
- contactPhysics->ks = jointShearStiffness*2.*R1*R2/(R1+R2);
- contactPhysics->dilationAngle = jointDilationAngle;
- contactPhysics->tanDilationAngle = std::tan(contactPhysics->dilationAngle);
+ Real jf1 = yade1->jointFrictionAngle;
+ Real jf2 = yade2->jointFrictionAngle;
+ Real jkn1 = yade1->jointNormalStiffness;
+ Real jkn2 = yade2->jointNormalStiffness;
+ Real jks1 = yade1->jointShearStiffness;
+ Real jks2 = yade2->jointShearStiffness;
+ Real jdil1 = yade1->jointDilationAngle;
+ Real jdil2 = yade2->jointDilationAngle;
+ Real jcoh1 = yade1->jointCohesion;
+ Real jcoh2 = yade2->jointCohesion;
+ Real jSigT1 = yade1->jointTensileStrength;
+ Real jSigT2 = yade2->jointTensileStrength;
+
+ //contactPhysics->frictionAngle = (jointFrictionAngle>=0 ? jointFrictionAngle : std::min(f1,f2));
+ contactPhysics->tanFrictionAngle = std::tan(std::min(jf1,jf2));
+
+ //contactPhysics->kn = jointNormalStiffness*2.*R1*R2/(R1+R2); // very first expression from Luc
+ //contactPhysics->kn = (jkn1+jkn2)/2.0*2.*R1*R2/(R1+R2); // after putting jointNormalStiffness in material
+ contactPhysics->kn = ( jkn1*Mathr::PI*pow(R1,2) + jkn2*Mathr::PI*pow(R2,2) )/2.0; // for a size independant expression
+
+ //contactPhysics->ks = jointShearStiffness*2.*R1*R2/(R1+R2);
+ //contactPhysics->ks = (jks1+jks2)/2.0*2.*R1*R2/(R1+R2);
+ contactPhysics->ks = ( jks1*Mathr::PI*pow(R1,2) + jks2*Mathr::PI*pow(R2,2) )/2.0; // for a size independant expression
+
+ contactPhysics->tanDilationAngle = std::tan(std::min(jdil1,jdil2));
///to set if the contact is cohesive or not
- if ((scene->iter < cohesiveTresholdIteration) && (jointCohesion>0 || jointTensileStrength>0)) {
+ if ((scene->iter < cohesiveTresholdIteration) && (std::min(jcoh1,jcoh2)>0 || std::min(jSigT1,jSigT2)/2.0>0)) {
contactPhysics->isCohesive=true;
- JCFpmState* st1=dynamic_cast<JCFpmState*>(Body::byId(interaction->getId1(),scene)->state.get());
st1->noIniLinks++;
- JCFpmState* st2=dynamic_cast<JCFpmState*>(Body::byId(interaction->getId2(),scene)->state.get());
st2->noIniLinks++;
}
else { contactPhysics->isCohesive=false; contactPhysics->FnMax=0; contactPhysics->FsMax=0; }
if ( contactPhysics->isCohesive ) {
- contactPhysics->FnMax = jointTensileStrength*contactPhysics->crossSection;
- contactPhysics->FsMax = jointCohesion*contactPhysics->crossSection;
+ contactPhysics->FnMax = std::min(jSigT1,jSigT2)*contactPhysics->crossSection;
+ contactPhysics->FsMax = std::min(jcoh1,jcoh2)*contactPhysics->crossSection;
}
}
interaction->phys = contactPhysics;
=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.hpp'
--- pkg/dem/JointedCohesiveFrictionalPM.hpp 2013-11-04 14:03:39 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.hpp 2013-11-04 14:20:24 +0000
@@ -9,12 +9,17 @@
/** This class holds information associated with each body state*/
class JCFpmState: public State {
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(JCFpmState,State,"JCFpm state information about each body.\n\nNone of that is used for computation (at least not now), only for post-processing.",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(JCFpmState,State,"JCFpm state information about each body.",
((int,tensBreak,0,,"number of tensile breakages. [-]"))
((int,shearBreak,0,,"number of shear breakages. [-]"))
((int,noIniLinks,0,,"number of initial cohesive interactions. [-]"))
- ((Real,tensBreakRel,0,,"relative number of tensile breakages (compared with noIniLinks). [-]"))
- ((Real,shearBreakRel,0,,"relative number of shear breakages (compared with noIniLinks). [-]"))
+ ((Real,tensBreakRel,0,,"relative number (in [0;1]) of tensile breakages (compared with noIniLinks). [-]"))
+ ((Real,shearBreakRel,0,,"relative number (in [0;1]) of shear breakages (compared with noIniLinks). [-]"))
+ ((bool,onJoint,false,,"identifies if the particle is on a joint surface."))
+ ((int,joint,0,,"indicates the number of joint surfaces to which the particle belongs (0-> no joint, 1->1 joint, etc..). [-]"))
+ ((Vector3r,jointNormal1,Vector3r::Zero(),,"specifies the normal direction to the joint plane 1. Rk: the ideal here would be to create a vector of vector wich size is defined by the joint integer (as much joint normals as joints). However, it needs to make the pushback function works with python since joint detection is done through a python script. lines 272 to 312 of cpp file should therefore be adapted. [-]"))
+ ((Vector3r,jointNormal2,Vector3r::Zero(),,"specifies the normal direction to the joint plane 2. [-]"))
+ ((Vector3r,jointNormal3,Vector3r::Zero(),,"specifies the normal direction to the joint plane 3. [-]"))
,
createIndex();
);
@@ -30,11 +35,16 @@
YADE_CLASS_BASE_DOC_ATTRS_CTOR(JCFpmMat,FrictMat,"possibly jointed cohesive frictional material, for use with other JCFpm classes",
((int,type,0,,"if particles of two different types interact, it will be with friction only (no cohesion).[-]"))
- ((bool,onJoint,false,,"identifies if the particle is on a joint surface. [-]"))
- ((int,joint,0,,"indicates the number of joint surfaces to which the particle belongs (0-> no joint, 1->1 joint, etc..). [-]"))
- ((Vector3r,jointNormal1,Vector3r::Zero(),,"specifies the normal direction to the joint plane 1. Rk: the ideal here would be to create a vector of vector wich size is defined by the joint integer (as much joint normals as joints). However, it needs to make the pushback function works with python since joint detection is done through a python script. lines 272 to 312 of cpp file should therefore be adapted. [-]"))
- ((Vector3r,jointNormal2,Vector3r::Zero(),,"specifies the normal direction to the joint plane 2. [-]"))
- ((Vector3r,jointNormal3,Vector3r::Zero(),,"specifies the normal direction to the joint plane 3. [-]")),
+ //((Real,alpha,0.,,"defines the shear to normal stiffness ratio ks/kn in the matrix."))
+ ((Real,tensileStrength,0.,,"defines the maximum admissible normal force in traction in the matrix (:yref:`JCFpmPhys.FnMax`=tensileStrength*:yref:`JCFpmPhys.crossSection`). [Pa]"))
+ ((Real,cohesion,0.,,"defines the maximum admissible tangential force in shear, for Fn=0, in the matrix (:yref:`JCFpmPhys.FsMax`=cohesion*:yref:`JCFpmPhys.crossSection`). [Pa]"))
+ ((Real,jointNormalStiffness,0.,,"defines the normal stiffness on the joint surface. [Pa/m]"))
+ ((Real,jointShearStiffness,0.,,"defines the shear stiffness on the joint surface. [Pa/m]"))
+ ((Real,jointTensileStrength,0.,,"defines the :yref:`maximum admissible normal force in traction<JCFpmPhys.FnMax>` on the joint surface. [Pa]"))
+ ((Real,jointCohesion,0.,,"defines the :yref:`maximum admissible tangential force in shear<JCFpmPhys.FsMax>`, for Fn=0, on the joint surface. [Pa]"))
+ ((Real,jointFrictionAngle,-1,,"defines Coulomb friction on the joint surface. [rad]"))
+ ((Real,jointDilationAngle,0,,"defines the dilatancy of the joint surface (only valid for :yref:`smooth contact logic<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.smoothJoint`). [rad]"))
+ ,
createIndex();
);
REGISTER_CLASS_INDEX(JCFpmMat,FrictMat);
@@ -48,18 +58,17 @@
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(JCFpmPhys,NormShearPhys,"representation of a single interaction of the JCFpm type, storage for relevant parameters",
((Real,initD,0,,"equilibrium distance for interacting particles. Computed as the interparticular distance at first contact detection."))
- ((bool,isCohesive,false,,"If false, particles interact in a frictional way. If true, particles are bonded regarding the given cohesion and tensileStrength. [-]"))
- ((bool,more,false,,"specifies if the interaction is crossed by more than 3 joints. If true, interaction is deleted (temporary solution). [-]"))
- ((bool,isOnJoint,false,,"if true, particles interact as part of a joint (see Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM for details). [-]"))
- ((Real,frictionAngle,0,,"defines Coulomb friction. [rad]"))
- ((Real,tanFrictionAngle,0,,"tangent of frictionAngle. [-]"))
+ ((bool,isCohesive,false,,"If false, particles interact in a frictional way. If true, particles are bonded regarding the given :yref:`cohesion<JCFpmMat.cohesion>` and :yref:`tensile strength<JCFpmMat.tensileStrength>` (or their jointed variants)."))
+ ((bool,more,false,,"specifies if the interaction is crossed by more than 3 joints. If true, interaction is deleted (temporary solution)."))
+ ((bool,isOnJoint,false,,"defined as true when both interacting particles are :yref:`on joint<JCFpmState.onJoint>` and are in opposite sides of the joint surface. In this case, mechanical parameters of the interaction are derived from the ''joint...'' material properties of the particles, and normal of the interaction is re-oriented (see also :yref:`Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM`)."))
+ ((Real,tanFrictionAngle,0,,"tangent of Coulomb friction angle for this interaction (auto. computed). [-]"))
((Real,crossSection,0,,"crossSection=pi*Rmin^2. [m2]"))
- ((Real,FnMax,0,,"defines the maximum admissible normal force in traction FnMax=tensileStrength*crossSection. [N]"))
- ((Real,FsMax,0,,"defines the maximum admissible tangential force in shear FsMax=cohesion*crossSection. [N]"))
- ((Vector3r,jointNormal,Vector3r::Zero(),,"Normal direction to the joint."))
+ ((Real,FnMax,0,,"computed from :yref:`tensile strength<JCFpmMat.tensileStrength>` (or joint variant) to define the maximum admissible normal force in traction. [N]"))
+ ((Real,FsMax,0,,"computed from :yref:`cohesion<JCFpmMat.cohesion>` (or jointCohesion) to define the maximum admissible tangential force in shear, for Fn=0. [N]"))
+ ((Vector3r,jointNormal,Vector3r::Zero(),,"normal direction to the joint, deduced from e.g. <JCFpmState.jointNormal1>."))
((Real,jointCumulativeSliding,0,,"sliding distance for particles interacting on a joint. Used, when :yref:`<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.smoothJoint>` is true, to take into account dilatancy due to shearing. [-]"))
- ((Real,dilationAngle,0,,"defines the dilatancy of the joint surface. [rad]"))
- ((Real,tanDilationAngle,0,,"tangent of dilationAngle. [-]"))
+ //((Real,dilationAngle,0,,"defines the dilatancy of the joint surface. [rad]"))
+ ((Real,tanDilationAngle,0,,"tangent of the angle defining the dilatancy of the joint surface (auto. computed from :yref:`JCFpmMat.jointDilationAngle`). [-]"))
((Real,dilation,0,,"defines the normal displacement in the joint after sliding treshold. [m]"))
,
createIndex();
@@ -80,15 +89,6 @@
YADE_CLASS_BASE_DOC_ATTRS(Ip2_JCFpmMat_JCFpmMat_JCFpmPhys,IPhysFunctor,"converts 2 JCFpmMat instances to JCFpmPhys with corresponding parameters.",
((int,cohesiveTresholdIteration,1,,"should new contacts be cohesive? They will before this iter, they won't afterward."))
- ((Real,alpha,0.,,"defines the shear to normal stiffness ratio ks/kn in the matrix."))
- ((Real,tensileStrength,0.,,"defines the maximum admissible normal force in traction FnMax=tensileStrength*crossSection in the matrix. [Pa]"))
- ((Real,cohesion,0.,,"defines the maximum admissible tangential force in shear FsMax=cohesion*crossSection in the matrix. [Pa]"))
- ((Real,jointNormalStiffness,0.,,"defines the normal stiffness on the joint surface. [Pa]"))
- ((Real,jointShearStiffness,0.,,"defines the shear stiffness on the joint surface. [Pa]"))
- ((Real,jointTensileStrength,0.,,"defines the maximum admissible normal force in traction FnMax=tensileStrength*crossSection on the joint surface. [Pa]"))
- ((Real,jointCohesion,0.,,"defines the maximum admissible tangential force in shear FsMax=cohesion*crossSection on the joint surface. [Pa]"))
- ((Real,jointFrictionAngle,-1,,"defines Coulomb friction on the joint surface. [rad]"))
- ((Real,jointDilationAngle,0,,"defines the dilatancy of the joint surface (only valid for smooth contact logic). [rad]"))
);
};
@@ -101,10 +101,10 @@
FUNCTOR2D(ScGeom,JCFpmPhys);
YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM,LawFunctor,"Interaction law for cohesive frictional material, e.g. rock, possibly presenting joint surfaces. Joint surfaces can be defined in a preprocessing phase through .stl meshes (see ref for details of the procedure), and can be mechanically described with a smooth contact logic [Ivars2011]_ (implemented in Yade in [Scholtes2012]_).",
- ((bool,smoothJoint,false,,"if true, particles belonging to joint surface have a smooth contact logic [Ivars2011]_, [Scholtes2012]_."))
- ((bool,recordCracks,false,,"if true a text file cracksKey.txt (see below) will be created (iteration, position, type (tensile or shear), cross section and contact normal)."))
- ((string,Key,"",,"string specifying the name of saved file 'cracks.....txt'"))
- ((bool,cracksFileExist,false,,"If true, text file already exists and its content won't be reset."))
+ ((bool,smoothJoint,false,,"if true, interactions of particles belonging to joint surface (:yref:`JCFpmPhys.isOnJoint`) are handled according to a smooth contact logic [Ivars2011]_, [Scholtes2012]_."))
+ ((bool,recordCracks,false,,"if true a text file cracksKey.txt (see below) will be created (apparition iteration, position, type (tensile or shear), cross section and contact normal)."))
+ ((string,Key,"",,"string specifying the name of saved file 'cracks___.txt', when :yref:`recordCracks<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.`recordCracks>` is true."))
+ ((bool,cracksFileExist,false,,"if true (and if :yref:`recordCracks<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.`recordCracks>`), data are appended to an existing 'cracks...' text file; otherwise its content is reset."))
);
DECLARE_LOGGER;
};
=== modified file 'py/export.py'
--- py/export.py 2013-10-07 14:56:14 +0000
+++ py/export.py 2013-11-04 14:20:24 +0000
@@ -47,7 +47,7 @@
elif (format=='id_x_y_z_r_matId'):
output+=('%d\t%g\t%g\t%g\t%g\t%d\n'%(b.id,b.state.pos[0],b.state.pos[1],b.state.pos[2],b.shape.radius,b.material.id))
elif (format=='jointedPM'):
- output+=('%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n'%(b.id,b.mat.onJoint,b.mat.joint,b.mat.jointNormal1[0],b.mat.jointNormal1[1],b.mat.jointNormal1[2],b.mat.jointNormal2[0],b.mat.jointNormal2[1],b.mat.jointNormal2[2],b.mat.jointNormal3[0],b.mat.jointNormal3[1],b.mat.jointNormal3[2]))
+ output+=('%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n'%(b.id,b.state.onJoint,b.state.joint,b.state.jointNormal1[0],b.state.jointNormal1[1],b.state.jointNormal1[2],b.state.jointNormal2[0],b.state.jointNormal2[1],b.state.jointNormal2[2],b.state.jointNormal3[0],b.state.jointNormal3[1],b.state.jointNormal3[2]))
elif (format=='liggghts_in'):
output+=('%g %g %g %g %g %g %g\n'%(count+1,b.mask,b.shape.radius,b.material.density,b.state.pos[0],b.state.pos[1],b.state.pos[2]))
outputVel+=('%g %g %g %g %g %g %g\n'%(count+1,b.state.vel[0],b.state.vel[1],b.state.vel[2],b.state.angVel[0],b.state.angVel[1],b.state.angVel[2]))
