New question #683679 on Yade:
https://answers.launchpad.net/yade/+question/683679

Hello:
I'm trying to generate a cylinder sample with target porosity, the MWE comes 
from [1], after I run the MWE, I click "show 3D" then click "Reference" or 
"Center" to show the model in different direction, there will be:
In [1]: Scene radius must be positive - Ignoring value
With this notification, the scene moves strange when I adjust it through mouse.
My yade version is 2018.02b with Ubuntu 18.04, here is the MWE:
################
from yade import pack
num_spheres=1000
key='_triax_base_'
targetPorosity = 0.45
compFricDegree = 30
finalFricDegree = 30
rate=-0.02
damp=0.2
stabilityThreshold=0.01
young=5e6
mn,mx=Vector3(0,0,0),Vector3(0.07,0.07,0.14)
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
Gl1_Sphere.quality=3

triax=TriaxialStressController(
        maxMultiplier=1.+2e4/young,
        finalMaxMultiplier=1.+2e3/young,
        thickness = 0,
        stressMask = 7,
        internalCompaction=False,
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        triax,
        TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+key),
        newton
]

triax.goal1=triax.goal2=triax.goal3=-10000

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001:
    break

print "###      Isotropic state saved      ###"

import sys
while triax.porosity>targetPorosity:
        compFricDegree = 0.95*compFricDegree
        setContactFriction(radians(compFricDegree))
        print "\r compFrictionDegree: ",compFricDegree," 
porosity:",triax.porosity,
        sys.stdout.flush()
        O.run(500,1)

print "###    Compacted state saved      ###"
#
aabb=utils.aabbExtrema()

predX=(aabb[0][0]+aabb[1][0])/2
predY=(aabb[0][1]+aabb[1][1])/2
predR=min((aabb[1][0]-aabb[0][0])/2,(aabb[1][1]-aabb[0][1])/2)

pred=pack.inCylinder((predX,predY,0),(predX,predY,0.14),predR)
# ### erase walls and spheres which is not in cylinder
def cylinderShape():
        for i in O.bodies:
                if not isinstance(i.shape,Sphere):
                        O.bodies.erase(i.id)
                if isinstance(i.shape,Sphere):
                                if not pred(i.state.pos, i.shape.radius):
                                        O.bodies.erase(i.id)
cylinderShape()
###############
Could you please help me with this problem?
Thanks in advance!
Leonard

[1]https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.

_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to     : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Reply via email to