# 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]
##########################################################################
##########################################################################
import itertools
from numpy import *
import numpy as np
import math
young=4e8
finalcompFricDegree=19.5
stabilityThreshold=0.02
############################################# sphere particles material properties ########################
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
for f in facets:
n = f.shape.normal
a = f.shape.area
##############################################################################################################
################################### 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
for f in facetstop:
n = f.shape.normal
a = f.shape.area
########################################################################################################
################################## 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
for f in facetsbottom:
n = f.shape.normal
a = f.shape.area
########################################################################################################
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=1000,command="checkunbalanced()"),
]
def checkunbalanced():
unb=unbalancedForce()
print 'unbalanced force: ',unb
if unb<stabilityThreshold:
O.pause()
export.text("cluster-final.txt")
O.saveTmp()
####################################################################################
####################################################################################

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
 Revision history for this message Klaus Thoeni (klaus.thoeni) said on 2021-02-18: #1

Hi,

facets in Yade are massless and non-dynamic. You cannot apply forces, you can only apply velicities/displacements. If you need dynamic facets consider using pFacets.

Klaus

 Revision history for this message Jan Stránský (honzik) said on 2021-02-18: #2

@ytang:

> the following code doesn't show any errors.

it does show this error:
IOError: [Errno 2] No such file or directory: 'cluster-adddd.txt'

> why the facets wall doesn't move?

The order of engines does matter.
In your case, the forces are reset before being applied (you set them by addForcestop, but then they are reset by ForceResetter).
Put addForcestop PyRunner **before** NewtonIntegrator (as it is in [2])

@Klaus
> facets in Yade are massless and non-dynamic. You cannot apply forces
this is the default state. By setting mass (and inertia) and unblocking desired DOFs, you can easily have dynamic "plain" facets.

But yes, pfacets may be a better choice for various reasons

cheers
Jan

 Revision history for this message ytang (ytang116) said on 2021-02-19: #3

Hi Jan, Klaus,

For Jan,

the file is in Dropbox.

You are right, I need to put the addForce function before the NewtonIntegrator.
Besides, I think it's very hard to control the facet, it's very to split it into pieces. So, I think I need to change to Pfacet.

best,
Yong

 Revision history for this message ytang (ytang116) said on 2021-02-19: #4

Thanks Jan Stránský, that solved my question.

 Revision history for this message Jan Stránský (honzik) said on 2021-02-20: #5

> the file is in Dropbox.

yes, it is clear from the link address :-) Nevertheless, it could not be publicly accessed (it is possible now), probably you will delete is some time etc. that's why there is "Please, no external links!" section in [4]

cheers
Jan