applying force on facet

Hi All,

I'm trying to simulate the CPT [1] in a deep location. In order to simulate the confining pressure state, I use the method mentioned here [2].

The basic idea is as follows:

(1) import the particles;
(2) create the top, bottom, and cylindrical walls by facets with mass
(3) give the facets with some forces.
In the end, I just want to test the unbalanced force. The result shows that the unbalanced force doesn't change at all even though I set the gravity force. Besides, it seems the facets don't move to reach the target pressure. (you can see these two situations here
[3] )
######################################################################
the following code doesn't show any errors.
My questions are:

(1) why the unbalanced force doesn't change? [3]
(2) why the facets wall doesn't move? [3]
##########################################################################
##########################################################################
from yade import ymport ,pack,export,geom,plot,export,utils
import itertools
from numpy import *
import numpy as np
import math
young=4e8
finalcompFricDegree=19.5
stabilityThreshold=0.02
############################################# sphere particles material properties ########################
sphereMat=O.materials.append(CohFrictMat(young=young,poisson=0.3,frictionAngle=radians(finalcompFricDegree),isCohesive=False,alphaKr=0.2,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True,density=2648))
O.bodies.append(ymport.text("cluster-adddd.txt",material = sphereMat)) #### the particles are here[3]
############################################################################################################
#################################
frictMat = O.materials.append(FrictMat(young=4e8,poisson=0.3,frictionAngle=19.5))
#######################################################################################################
############################################### cylindrical facets ####################################
width = 0.4 ; height = 0.5 ; preStress = -3e5
#########################facets division
nw = 25 ; nh = 15
############################################### cylindrical facets ####################
facets = []
rCyl2 = .5*width / cos(pi/float(nw))
for r in range(nw):
 for h in range(nh):
  v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw))+0.2, rCyl2*sin(2*pi*(r+0)/float(nw))+0.2, height*(h+0)/float(nh) )
  v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw))+0.2, rCyl2*sin(2*pi*(r+1)/float(nw))+0.2, height*(h+0)/float(nh) )
  v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw))+0.2, rCyl2*sin(2*pi*(r+1)/float(nw))+0.2, height*(h+1)/float(nh) )
  v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw))+0.2, rCyl2*sin(2*pi*(r+0)/float(nw))+0.2, height*(h+1)/float(nh) )
  f1 = facet((v1,v2,v3),color=(0,0,1),material=frictMat)
  f2 = facet((v1,v3,v4),color=(0,0,1),material=frictMat)
  facets.extend((f1,f2))
O.bodies.append(facets)
mass = O.bodies[0].state.mass
for f in facets:
 f.state.mass = mass
 f.state.blockedDOFs = 'XYZz'
############################################# apply prestress to facets
def addForces():
 for f in facets:
  n = f.shape.normal
  a = f.shape.area
  O.forces.addF(f.id,preStress*a*n)
##############################################################################################################
################################### top facets ##############################################################
widthtop = 0.5 ; heighttop = 0.45 ; preStresstop = -3e5
######################################################## facets division
nxtop = 25 ; nytop = 25
################################################
facetstop = []
for nx in range(nxtop):
 for ny in range(nytop):
  v11 = Vector3( (nx+0)*widthtop/nxtop-0.05, (ny+0)*widthtop/nytop-0.05, heighttop )
  v12 = Vector3( (nx+1)*widthtop/nxtop-0.05, (ny+0)*widthtop/nytop-0.05, heighttop )
  v13 = Vector3( (nx+1)*widthtop/nxtop-0.05, (ny+1)*widthtop/nytop-0.05, heighttop )
  v14 = Vector3( (nx+0)*widthtop/nxtop-0.05, (ny+1)*widthtop/nytop-0.05, heighttop )
  f11 = facet((v11,v12,v13),color=(0,0,1),material=frictMat)
  f12 = facet((v11,v13,v14),color=(0,0,1),material=frictMat)
  facetstop.extend((f11,f12))
O.bodies.append(facetstop)
masstop = O.bodies[0].state.mass
for f in facetstop:
 f.state.mass = masstop
 f.state.blockedDOFs = 'XYZxy'
############################################# apply prestress to facets
def addForcestop():
 for f in facetstop:
  n = f.shape.normal
  a = f.shape.area
  O.forces.addF(f.id,preStresstop*a*n)
########################################################################################################
################################## bottom facets########################################################
widthbottom = 0.5 ; heightbottom = 0.0 ; preStressbottom = 3e5
#########################facets division
nxbottom = 25 ; nybottom = 25
facetsbottom = []
for nx in range(nxbottom):
 for ny in range(nybottom):
  v11 = Vector3( (nx+0)*widthbottom/nxbottom-0.05, (ny+0)*widthbottom/nybottom-0.05, heightbottom )
  v12 = Vector3( (nx+1)*widthbottom/nxbottom-0.05, (ny+0)*widthbottom/nybottom-0.05, heightbottom )
  v13 = Vector3( (nx+1)*widthbottom/nxbottom-0.05, (ny+1)*widthbottom/nybottom-0.05, heightbottom )
  v14 = Vector3( (nx+0)*widthbottom/nxbottom-0.05, (ny+1)*widthbottom/nybottom-0.05, heightbottom )
  f11 = facet((v11,v12,v13),color=(0,0,1),material=frictMat)
  f12 = facet((v11,v13,v14),color=(0,0,1),material=frictMat)
  facetsbottom.extend((f11,f12))
O.bodies.append(facetsbottom)
massbottom = O.bodies[0].state.mass
for f in facetsbottom:
 f.state.mass = massbottom
 f.state.blockedDOFs = 'XYZxy'
############################################# apply prestress to facets
def addForcesbottom():
 for f in facetsbottom:
  n = f.shape.normal
  a = f.shape.area
  O.forces.addF(f.id,preStressbottom*a*n)
########################################################################################################
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom()],
   [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
     [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True,label='cohesiveLaw')]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.0),
 VTKRecorder(fileName='3d-vtk-',recorders=['all','bstresses'],iterPeriod=20000),
 ############################################################################
 ############################################################################
 PyRunner(iterPeriod=1,command="addForces()"),
 PyRunner(iterPeriod=1,command="addForcestop()"),
 PyRunner(iterPeriod=1,command="addForcesbottom()"),
 PyRunner(iterPeriod=1000,command="checkunbalanced()"),
]
def checkunbalanced():
 unb=unbalancedForce()
 print 'unbalanced force: ',unb
 if unb<stabilityThreshold:
  O.pause()
export.text("cluster-final.txt")
O.save('make_sample_cluster_final-confining.yade.gz' )
O.saveTmp()
####################################################################################
####################################################################################

references:
[1] https://en.wikipedia.org/wiki/Cone_penetration_test
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/concrete/triax.py
[3] https://www.dropbox.com/sh/r8j8cpviz742y21/AAAMJ9opuVFWkiPLyioeetJ_a?dl=0