[YADE] Cube sample generation

  # -*- coding: utf-8 -*-
from yade import pack
from yade import utils
import math
import os
#from utils import *

num_spheres= 10000
mn,mx=Vector3(-0.1,-0.1,-0.1),Vector3(0.1,0.1,0.1) # Top a bottom point of the cube
thick = 0.01
compFricDegree = table.compFricDegree
finalFricDegree=35
#targetPorosity=0.382
rate=0.01
damp=0.06
stabilityThreshold=1e-3 # leave it too low can have a super stable state, but will eat onion for waiting so long

key='_sample' # just used for adding name to the output, kind of boring

## create material #0, which will be used as default
O.materials.append(FrictMat(young=356e6,poisson=0.42,frictionAngle=radians(compFricDegree),density=3000,label='spheres')) # hat
O.materials.append(FrictMat(young=356e6,poisson=0.5,frictionAngle=0,density=0,label='walls')) # boundary

## create walls around the packing
walls=utils.aabbWalls([mn,mx],thickness=thick,oversizeFactor=3,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
#psdSizes=[0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.0095] # (sizes or radii of the grains vary from 2mm to 9.5mm)
#psdCumm=[0.01,0.09,0.25,0.50,0.69,0.90,0.95,1] # the correspondent amount (percentage) of each diameter
#psdSizes,psdCumm=[0.001978,0.00218,0.00249,0.002744,0.002894,0.002999,0.003162,0.003305,0.003455,0.00358,0.003709,0.003911,0.004016,0.004273,0.004427,0.004669,0.004968,0.005193,0.005476,0.005826,0.006036,0.00631,0.006595,0.006833,0.007079,0.007335,0.007532,0.007943,0.008526,0.0092],[0.001,0.002,0.005,0.026385,0.047493,0.068602,0.092348,0.12372,0.150396,0.174142,0.200528,0.242744,0.261214,0.311346,0.356201,0.401055,0.469657,0.522427,0.583113,0.654354,0.71504,0.770449,0.823219,0.878628,0.926121,0.968338,0.98153,0.992084,0.997361,1]
#---------------------------------------------
#sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,psdSizes,psdCumm,False,seed=1)
if num_spheres==1000:
  sp.makeCloud(mn,mx,0.005,0.65,num_spheres,False,0.8,seed=1) #initial 0.001 for 1000 particle use 0.005
elif num_spheres==10000:
  sp.makeCloud(mn,mx,0.001,0.65,num_spheres,False,0.8,seed=1)
  #sp.makeCloud(mn,mx,-1,0,-1,False, 0.95,psdSizes,psdCumm,False,seed=1)
sp.toSimulation(material='spheres')

volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
mean_rad = pow(0.09*volume/num_spheres,0.3333)

#clumps=False
#if clumps:
#   c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
#   sp.makeClumpCloud((-0.24,0.24,-0.24),(0.24,0.24,0.24),[c1],periodic=False)
#   O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
#   standalone,clumps=sp.getClumps()
#   for clump in clumps:
#       O.bodies.clump(clump)
#       for i in clump[1:]: O.bodies[i].shape.color=O.bodies[clump[0]].shape.color
#   #sp.toSimulation()
#else:
#   O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

O.dt=.5*utils.PWaveTimeStep() # initial timestep, to not explode right away
O.usesTimeStepper=True

triax=TriaxialStressController(
    maxMultiplier=1.001,
    finalMaxMultiplier=1.0001,
    thickness = thick,
    radiusControlInterval=10
#   stressMask = 7,
#   max_vel=0.01,
#   strainDamping=0.99, #0.99
#   stressDamping=0.25
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],verletDist=-mean_rad*0.06),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack(label='law')]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
    triax,
    TriaxialStateRecorder(iterPeriod=50,file='WallStresses'+key,addIterNum=False,initRun=False),
    newton
] #timestep update interval is 100 by default

#O.save('initial'+key+'.xml')
#Display spheres with 2 colors for seeing rotations better
#Gl1_Sphere.stripes=0
#if nRead==0: yade.qt.Controller(), yade.qt.View()
print 'Number of elements: ', len(O.bodies)
print 'Box Volume: ',  triax.boxVolume
print 'Box Volume calculated: ', volume
#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=1
#yade.qt.Controller(), yade.qt.View()

#########################################
# Phase 1: ISOTROPIC GENERATOR OF 20kPa #
#########################################
while 1:
  triax.internalCompaction=True # Growing the particles is what we use in this step
  setContactFriction(radians(compFricDegree))
  triax.stressmask=7
  triax.goal1=10000
  triax.goal2=10000
  triax.goal3=10000
  O.run(10,True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
  print 'unbalanced force:',unb,' mean stress: ',meanS
#  print 'void ratio=',triax.porosity/(1-triax.porosity), 'porosity=', triax.porosity
  print 'mean stress engine', triax.meanStress, 'kineticE=', utils.kineticEnergy()
  print '-----State_01: Isotropic compression 10kPa-----'
  if unb<stabilityThreshold and abs(meanS-table.isoForce)/table.isoForce<0.001:
    break

O.save('confinedState20'+key+'.xml') # remember this impotant part, will load is later
print "##   Isotropic state saved (20 kPa)  ##"
print 'current porosity=',triax.porosity
print 'current void ratio=',triax.porosity/(1-triax.porosity)