TriaxialStressController and Capillary

from yade import pack, plot
from yade import utils
import math
import os

num_spheres=500 # 10000 spheres
mn,mx=Vector3(-0.1,-0.1,-0.1),Vector3(0.1,0.1,0.1)
thick = 0 # thickness used for walls (not so important)
compFricDegree = 1 # to generate a super dense specimen
wall_velocity = 0.5

rate=0.01
damp=0.06

local_damping  = 0.006

## create material
O.materials.append(FrictMat(young=50e6,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=50e6,poisson=0.5,frictionAngle=0,density=0,label='walls'))

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

sp=pack.SpherePack()
sp.makeCloud(mn,mx,0.005,0.3,num_spheres,False,0.8,seed=1)
O.bodies.append([sphere(c,r,material='spheres') for c,r in sp])

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

O.dt=.5*utils.PWaveTimeStep()
O.usesTimeStepper=True

triax=TriaxialStressController(
    maxMultiplier         =1.001,
    finalMaxMultiplier    =1.000001,
    thickness             = thick,
    radiusControlInterval =10,
    stressMask            = 7,
    max_vel               =0.01,
    strainDamping         =0.75,
    stressDamping         =0.01
)

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
            [Ip2_FrictMat_FrictMat_CapillaryPhys(label='capillary_phys')],
            [Law2_ScGeom_FrictPhys_CundallStrack(label='cundall_law')]  #for linear model only
        ),
        Law2_ScGeom_CapillaryPhys_Capillarity(capillaryPressure=1000.0,label='cap_law'),#for linear model only
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
        triax,
        NewtonIntegrator(damping=local_damping,gravity=(0,0,0))
    ]

cap_law.dead=False
cundall_law.neverErase=True
triax.internalCompaction = True

from yade import qt

triax.stressMask = 7
triax.goal1 = 5000
triax.goal2 = 5000
triax.goal3 = 5000
qt.View()