addForce didnt work for facet

Asked by onyourself on 2020-10-14

hi,

1) i generate a 2D triax specimen(x-z plane) of which the boundary is consisted with facet, then i add force on them with the code below, it just doesnt work.

2) whats the normal direction of a facet? eg : normal direction of a facet on x-y palne is (0,0,1) or (0,0,-1)? does this has something to do with the way how vertices were connected to a facet?

#########################
# facets division
nw = 240,
nh = 150,

#meterials########
O.materials.append(FrictMat(young=4.0e6,poisson=.3,frictionAngle=20,density=910
e+6,label='sphcylMat'))
O.materials.append(FrictMat(young=4.0e6,poisson=.3,frictionAngle=26,density=910
e+6,label='cylinderMat'))
O.materials.append(FrictMat(young=4.0e6,poisson=.15,frictionAngle=30,density=26
30e+6,label='sphereMat'))
O.materials.append(FrictMat(young=4.0e6,poisson=.15,frictionAngle=0,density=910
e+6,label='walls'))
#O.materials.append(FrictMat(young=4.0e6,poisson=.3,frictionAngle=26,density=91
0e+6,label='frictMat'))
frictMat = O.materials.append(FrictMat(
young=4.0e6,poisson=.3,frictionAngle=26,density=910e+6
))

facets = []
nw2 = nw/40
nw2 = 240
for r in xrange(nw2):
for h in xrange(nh):
v11 = Vector3( -.5*width + (r+0)*width/nw2, 5e-3, 0 )
v12 = Vector3( -.5*width + (r+1)*width/nw2, -5e-3, 0 )
v13 = Vector3( -.5*width + (r+1)*width/nw2, 5e-3, 0 )
v14 = Vector3( -.5*width + (r+0)*width/nw2, -5e-3, 0 )
f11 = facet((v11,v12,v13),color=(0,0,1),material=frictMat)
f12 = facet((v11,v13,v14),color=(0,0,1),material=frictMat)
v21 = Vector3( +.5*width, -5e-3, height*(h+0)/float(nh) )
v22 = Vector3( +.5*width, 5e-3, height*(h+0)/float(nh) )
v23 = Vector3( +.5*width, 5e-3, height*(h+1)/float(nh) )
v24 = Vector3( +.5*width, -5e-3, height*(h+1)/float(nh) )
f21 = facet((v21,v22,v23),color=(0,0,1),material=frictMat)
f22 = facet((v21,v23,v24),color=(0,0,1),material=frictMat)
v31 = Vector3( +.5*width - (r+0)*width/nw2, +5e-3, 250e-3 )
v32 = Vector3( +.5*width - (r+1)*width/nw2, +5e-3, 250e-3 )
v33 = Vector3( +.5*width - (r+1)*width/nw2, -5e-3, 250e-3 )
 v34 = Vector3( +.5*width - (r+0)*width/nw2, -5e-3, 250e-3 )
f31 = facet((v31,v32,v33),color=(0,0,1),material=frictMat)
f32 = facet((v31,v33,v34),color=(0,0,1),material=frictMat)
v41 = Vector3( -.5*width, -5e-3, height*(h+0)/float(nh) )
v42 = Vector3( -.5*width, 5e-3, height*(h+0)/float(nh) )
v43 = Vector3( -.5*width, 5e-3, height*(h+1)/float(nh) )
v44 = Vector3( -.5*width, -5e-3, height*(h+1)/float(nh) )
f41 = facet((v41,v42,v43),color=(0,0,1),material=frictMat)
f42 = facet((v41,v43,v44),color=(0,0,1),material=frictMat)
facets.extend((f11,f12,f21,f22,f31,f32,f41,f42))
O.bodies.append(facets)
mass = O.bodies[0].state.mass
for f in facets:
f.state.mass = mass
f.shape.color=(200,255,240)
f.state.blockedDOFs = 'XYZz'

def addForces():
 for f in facets:
  n = f.shape.normal
  a = f.shape.area
  O.forces.addF(f.id,preStress*a*n)
##########################
my whole code is a bit long so i upload related parts of it, thanks.

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-10-21
Last query:
2020-10-21
Last reply:
2020-10-16
Jan Stránský (honzik) said : #1

Hello,

1) please provide a MWE [1], W=working, reproducing your problem (but still keeping it M=minimal). Here, addForces is never called, so you do not apply the force actually and that is the reason why "it just doesnt work"...
There can be many reasons, hard even to guess without actual code:
- addForces / O.forces.addF is never used / is misused
- preStress is 0 (not defined in the provided code)
- it just works, just you wrongly interpret it (typically it is not "visible in GUI" but observable according to numerical results)
- ...

2) yes [2], in python it would be
normal = (v1-v0).cross(v2-v1).normalized() # vertices are v0,v1,v2
so the normal is "right-hand-ruled by v0,v1,v2"

cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/common/Facet.cpp#L36

onyourself (xxxe) said : #2

hi,

1) - it just works, just you wrongly interpret it (typically it is not "visible in GUI" but observable according to numerical results)

if so, maybe i just need more time to see the difference? i didnt print anything so the only way i can be sure that my code works is to check 3D view by my eyes…i thought it should be visible as my specimen is actual incompact.

2)got it, thanks, below is my code.

####################################
from yade.gridpfacet import *
from yade import pack, plot
from random import random
import numpy as np
from numpy import *
import math
import os

#parameters
rParticle = 0.14e-3#diameter0.28
rRelFuzz=.0005
mi,ma = (-500e-3,0.,0.),(500e-3,0,250e-3)
#nCyls,nSphs = 8000,28216
frictionAngleSph=30
frictionAngleCyl=30
width = 1000.e-3,
length = 10e-3,
height = 250.e-3,

# default parameters or from table
readParamsFromTable(noTableOk=True,
 # type of test ['cyl','cube']
 testType = 'cube',

 # material parameters
 young = 20e9,
 poisson = .2,
 frictionAngle = 1.2,
 sigmaT = 1.5e6,
 epsCrackOnset = 1e-4,
 relDuctility = 30,

 # prestress
 preStress = -3e9,
 # axial strain rate
 strainRate = -100,

 # assamlby parameters
 rParticle = .075e-3, #
 width = 1000e-3,
 height = 250e-3,
 bcCoeff = 5,

 # facets division
 nw = 240,
 nh = 150,

 # output specifications
 fileName = 'test',
 exportDir = '/tmp',
 runGnuplot = False,
 runInGui = True,
)
from yade.params.table import *

#meterials########
O.materials.append(FrictMat(young=4.0e6,poisson=.3,frictionAngle=20,density=910e+6,label='sphcylMat'))
O.materials.append(FrictMat(young=4.0e6,poisson=.3,frictionAngle=26,density=910e+6,label='cylinderMat'))
O.materials.append(FrictMat(young=4.0e6,poisson=.15,frictionAngle=30,density=2630e+6,label='sphereMat'))
O.materials.append(FrictMat(young=4.0e6,poisson=.15,frictionAngle=0,density=910e+6,label='walls'))
#O.materials.append(FrictMat(young=4.0e6,poisson=.3,frictionAngle=26,density=910e+6,label='frictMat'))
frictMat = O.materials.append(FrictMat(
 young=4.0e6,poisson=.3,frictionAngle=26,density=910e+6
))

##sphere
sp = yade.pack.SpherePack()
sp.makeCloud(mi,ma,porosity=0.7/1.7,psdSizes=[.00014,0.00016,0.00022,.0003,.00035],psdCumm=[0.,0.1,0.3,0.6,1.],num=35000)
#psdSizes=[0.01,0.0125,0.015,0.018,0.02,0.023,0.025,0.027,0.029,0.03],psdCumm=[0,0.05,0.15,0.38,0.54,0.77,0.9,0.97,0.99,1]) # makeCloud for spheres
spheres=sp.toSimulation(color=(0,0.5,0.7),material='sphereMat')
print "cheers"

num=0
vol=0
mass=0
area=0
for b in O.bodies:
 if isinstance(b.shape,Sphere):
     r=b.shape.radius
     x=b.state.pos[0]
     y=b.state.pos[1]
     z=b.state.pos[2]
     ID=b.id
     num+=1
     #vol+=4./3.*np.pi*r**3
     mass+=2.63*4./3.*np.pi*r**3
   area+=np.pi*r**2
print 'num,mass,voidRatio =',num,mass,((1000e-3*250e-3)-area)/area

length=6e-3
rfiber=0.023e-3/2
fibre_mass=np.pi*(rfiber)**2*length*0.91
Xw=0.5*1e-2 # fiber content
numFiber=Xw*mass/fibre_mass
print "the numer of fibres =",int(numFiber)

Ne=int(length/(0.28e-3/2))
print 'Ne=',Ne
nodesIds=[]
cylIds=[]
numFibre=0
#target_num=int(numFiber)
target_num=1200
deleted_sphere=0

####fiberfilewrite
fiberwr=open('fibers.txt','w')
cylNodes=[]

for n in range(2000):
 if numFibre<target_num:

  random.seed()
  l_x=random.uniform(-500e-3,500e-3)
  l_z=random.uniform(0,250e-3)
  x0=l_x
  z0=l_z

  u1 = random.random()
  u2 = random.random()*pi
  u3 = random.random()*pi
  a = sqrt(1-u1)
  b = sqrt(u1)
  t=((a*cos(u2))**2+(b*sin(u3))**2)**0.5
  cosphi=a*cos(u2)/t
  sinphi=b*sin(u3)/t

  hx=length*cosphi
  hz=length*sinphi

  x1=x0+hx
  y1=0
  z1=z0+hz

  if -500e-3<x1<500e-3 and 0<z1<250e-3 and -500e-3<x0<500e-3 and 0<z0<250e-3:
   if len(cylNodes)<=1:
    fiberwr.write(str(numFibre)+'\t'+str(x0)+'\t'+str(0)+'\t'+str(z0)+'\t'+str(x1)+'\t'+str(0)+'\t'+str(z1)+'\n')
    cylNodes.append([x0,z0,x1,z1])
    numFibre+=1
    vertices=[]
    for i in range(0, Ne+1):
     px=float(i)*hx/float(Ne)+x0; py=0; pz=float(i)*hz/float(Ne)+z0;
     vertices.append([px,py,pz])
    cylinderConnection(vertices,0.023e-3/2,nodesIds,cylIds,color=[0.5,0.5,0],highlight=False,intMaterial='cylinderMat',extMaterial='fMat')
    print "+1"
   else:
    valid=True
    #print len(cylNodes)
    for c in cylNodes:
     if(max(c[0],c[2])>=min(x0,x1)
     and max(x0,x1)>=min(c[0],c[2])
     and max(c[1],c[3])>=min(z0,z1)
     and max(z0,z1)>=min(c[1],c[3])):
      x12=(x1-x0)
      z12=(z1-z0)
      x13=(c[0]-x0)
      z13=(c[1]-z0)
      x14=(c[2]-x0)
      z14=(c[3]-z0)
      cross1=x12*z13-x13*z12
      cross2=x12*z14-x14*z12
      x34=c[2]-c[0]
      z34=c[3]-c[1]
      #x23=c[0]-x1
      #z23=c[1]-z1
      x32=x1-c[0]
      z32=z1-c[1]
      #x24=c[2]-x1
      #z24=c[3]-z1
      x31=x0-c[0]
      z31=z0-c[1]
      cross3=x34*z32-x32*z34
      cross4=x34*z31-x31*z34
      if(cross1*cross2<=0
      and cross3*cross4<=0):
       valid=False
       print "next"
       break
    if valid:
     fiberwr.write(str(numFibre)+'\t'+str(x0)+'\t'+str(0)+'\t'+str(z0)+'\t'+str(x1)+'\t'+str(0)+'\t'+str(z1)+'\n')
     cylNodes.append([x0,z0,x1,z1])
     numFibre+=1
     #print "++1"
     vertices=[]
     for i in range(0, Ne+1):
      px=float(i)*hx/float(Ne)+x0; py=0; pz=float(i)*hz/float(Ne)+z0;
      vertices.append([px,py,pz])
     cylinderConnection(vertices,0.023e-3/2,nodesIds,cylIds,color=[0.5,0.5,0],highlight=False,intMaterial='cylinderMat',extMaterial='fMat')

   ###spherefileWrite
   ##################################################################
   sphwr=open('spheres.txt','w')
   for b in O.bodies:
    if isinstance(b.shape,Sphere):
         r=b.shape.radius
     if r>0.012e-3:
          x=b.state.pos[0]
         y=b.state.pos[1]
          z=b.state.pos[2]
      ID=b.id
      k=(z1-z0)/(x1-x0)
      #kx-z+z0-kx0=0
      A=k
      B=-1
      C=z0-k*x0
      #if (A*x+B*z+C)/(A^2+B^2)^0.5<=(0.14e-3+0.023e-3/2):
      if abs(k*x-1*z+C)/(k**2+1)**0.5<=(r+0.0115e-3):
       O.bodies.erase(b.id)
       deleted_sphere+=1
      #else:
      sphwr.write(str(ID)+'\t'+str(x)+'\t'+str(y)+'\t'+str(z)+'\t'+str(r)+'\n')
   sphwr.close()
   ##################################################################
fiberwr.close()
print 'numFibre , deleted_sphere =',numFibre, deleted_sphere

for i in O.bodies:
        if isinstance(i.shape,Sphere):
         i.state.blockedDOFs="yXZ"
         #i.shape.color=(0.,0.5,0.7)

facets = []
nw2 = nw/40
nw2 = 240
for r in xrange(nw2):
 for h in xrange(nh):
  v11 = Vector3( -.5*width + (r+0)*width/nw2, 5e-3, 0 )
  v12 = Vector3( -.5*width + (r+1)*width/nw2, -5e-3, 0 )
  v13 = Vector3( -.5*width + (r+1)*width/nw2, 5e-3, 0 )
  v14 = Vector3( -.5*width + (r+0)*width/nw2, -5e-3, 0 )
  f11 = facet((v11,v12,v13),color=(0,0,1),material=frictMat)
  f12 = facet((v11,v13,v14),color=(0,0,1),material=frictMat)
  v21 = Vector3( +.5*width, -5e-3, height*(h+0)/float(nh) )
  v22 = Vector3( +.5*width, 5e-3, height*(h+0)/float(nh) )
  v23 = Vector3( +.5*width, 5e-3, height*(h+1)/float(nh) )
  v24 = Vector3( +.5*width, -5e-3, height*(h+1)/float(nh) )
  f21 = facet((v21,v22,v23),color=(0,0,1),material=frictMat)
  f22 = facet((v21,v23,v24),color=(0,0,1),material=frictMat)
  v31 = Vector3( +.5*width - (r+0)*width/nw2, +5e-3, 250e-3 )
  v32 = Vector3( +.5*width - (r+1)*width/nw2, +5e-3, 250e-3 )
  v33 = Vector3( +.5*width - (r+1)*width/nw2, -5e-3, 250e-3 )
  v34 = Vector3( +.5*width - (r+0)*width/nw2, -5e-3, 250e-3 )
  f31 = facet((v31,v32,v33),color=(0,0,1),material=frictMat)
  f32 = facet((v31,v33,v34),color=(0,0,1),material=frictMat)
  v41 = Vector3( -.5*width, -5e-3, height*(h+0)/float(nh) )
  v42 = Vector3( -.5*width, 5e-3, height*(h+0)/float(nh) )
  v43 = Vector3( -.5*width, 5e-3, height*(h+1)/float(nh) )
  v44 = Vector3( -.5*width, -5e-3, height*(h+1)/float(nh) )
  f41 = facet((v41,v42,v43),color=(0,0,1),material=frictMat)
  f42 = facet((v41,v43,v44),color=(0,0,1),material=frictMat)
  facets.extend((f11,f12,f21,f22,f31,f32,f41,f42))
O.bodies.append(facets)
mass = O.bodies[0].state.mass
for f in facets:
 f.state.mass = mass
 f.shape.color=(200,255,240)
 f.state.blockedDOFs = 'XYZz'
'''
#preStress
# 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)
'''
for f in facets:
 n = f.shape.normal
 a = f.shape.area
 #O.forces.addF(f.id,preStress*a*n)
 O.forces.setPermF(f.id,preStress*a*n)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  #Bo1_Box_Aabb(),
  Bo1_Sphere_Aabb(),
  Bo1_GridConnection_Aabb(),
  Bo1_Facet_Aabb(),
 ]),
 InteractionLoop([
  Ig2_Sphere_Sphere_ScGeom(),
  #Ig2_Box_Sphere_ScGeom(),
  Ig2_GridNode_GridNode_GridNodeGeom6D(),
  Ig2_Sphere_GridConnection_ScGridCoGeom(),
  Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
  Ig2_Facet_Sphere_ScGeom(),
 ],
 [
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False), # internal cylinder physics
  Ip2_FrictMat_FrictMat_FrictPhys() # physics for external interactions, i.e., cylinder-cylinder, sphere-sphere, cylinder-sphere
  #Ip2_FrictMat_CpmMat_FrictPhys(),
 ],
 [
  Law2_ScGeom_FrictPhys_CundallStrack(), # contact law for sphere-sphere
  Law2_ScGridCoGeom_FrictPhys_CundallStrack(), # contact law for cylinder-sphere
  Law2_ScGeom6D_CohFrictPhys_CohesionMoment(), # contact law for "internal" cylinder forces
  Law2_GridCoGridCoGeom_FrictPhys_CundallStrack() # contact law for cylinder-cylinder interaction
  #Law2_ScGeom_CpmPhys_Cpm(),
 ]
 ),
 GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.1,label='ts'),
 #triax,
 #TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.1,label='newton'),
 #PyRunner(iterPeriod=200,command='targetSpecimen()',label='target'),
 #PyRunner(iterPeriod=20,command='history()',label='recorder'),
]

######################################################

Jan Stránský (honzik) said : #3

Hello,

> if so, maybe i just need more time to see the difference? i didnt print anything so the only way i can be sure that my code works is to check 3D view by my eyes…i thought it should be visible as my specimen is actual incompact.

depneding on "to see" is a bad approach. The only way how to be sure is NOT checking 3D view, the only way how to be sure are the numbers in the model.

Thanks for the code, could you try to make it minimal and also not-so-long-running? I have started it a few minutes ago, still waiting for the initial configuration to finish..

O.forces.setPermF(f.id,preStress*a*n) is the right approach

> f.state.blockedDOFs = 'XYZz'

check if normal is not parallel to z, which is blocked

Otherwise, do
###
... your code

f = facets[0] # or any other valid index
print(f.state.pos)
O.step()
print(f.state.pos) # should be different from the first one
###

cheers
Jan

onyourself (xxxe) said : #4

hi,

> depneding on "to see" is a bad approach. The only way how to be sure is NOT checking 3D view, the only way how to be sure are the numbers in the model.

got it. then i change my code as you suggest(print(f.state.pos)). the positions printed twice are the same. please help me figure out why. i try my best to minimize my code.

############################################
from yade.gridpfacet import *
from yade import pack, plot
from random import random
from numpy import *
import math
import os

#parameters
rParticle = 0.14e-3#diameter0.28
rRelFuzz=.0005
mi,ma = (-500e-3,0.,0.),(500e-3,0,250e-3)
length = 10e-3,
width = int(1000e-3),
height = int(250e-3),

# prestress
preStress = -3e9,
# facets division
nw = int(240)
nh = int(150)

#meterials########
frictMat = O.materials.append(FrictMat(
 young=4.0e6,poisson=.3,frictionAngle=26,density=910e+6
))

facets = []
nw2 = (250)
nh = (50)
width = 1000e-3
height = 250e-3
for r in xrange(nw2):
 for h in xrange(nh):
                v11 = Vector3( -.5*width + (r+0)*width/nw2, 5e-3, 0 )
  v12 = Vector3( -.5*width + (r+1)*width/nw2, 5e-3, 0 )
  v13 = Vector3( -.5*width + (r+1)*width/nw2, -5e-3, 0 )
  v14 = Vector3( -.5*width + (r+0)*width/nw2, -5e-3, 0 )
  #print (v11, v12, v13, v14)
  f11 = facet((v11,v12,v13),color=(0,0,1),material=frictMat)
  f12 = facet((v11,v13,v14),color=(0,0,1),material=frictMat)
  v21 = Vector3( +.5*width, -5e-3, height*(h+0)/float(nh) )
  v22 = Vector3( +.5*width, 5e-3, height*(h+0)/float(nh) )
  v23 = Vector3( +.5*width, 5e-3, height*(h+1)/float(nh) )
  v24 = Vector3( +.5*width, -5e-3, height*(h+1)/float(nh) )
  #print (v21, v22, v23, v24)
  f21 = facet((v21,v22,v23),color=(0,0,1),material=frictMat)
  f22 = facet((v21,v23,v24),color=(0,0,1),material=frictMat)
  v31 = Vector3( +.5*width - (r+0)*width/nw2, +5e-3, 250e-3 )
  v32 = Vector3( +.5*width - (r+1)*width/nw2, +5e-3, 250e-3 )
  v33 = Vector3( +.5*width - (r+1)*width/nw2, -5e-3, 250e-3 )
  v34 = Vector3( +.5*width - (r+0)*width/nw2, -5e-3, 250e-3 )
  #print (v31, v32, v33, v34)
  f31 = facet((v31,v32,v33),color=(0,0,1),material=frictMat)
  f32 = facet((v31,v33,v34),color=(0,0,1),material=frictMat)
  v41 = Vector3( -.5*width, -5e-3, height*(h+0)/float(nh) )
  v42 = Vector3( -.5*width, 5e-3, height*(h+0)/float(nh) )
  v43 = Vector3( -.5*width, 5e-3, height*(h+1)/float(nh) )
  v44 = Vector3( -.5*width, -5e-3, height*(h+1)/float(nh) )
  #print (v41, v42, v43, v44)
  f41 = facet((v41,v42,v43),color=(0,0,1),material=frictMat)
  f42 = facet((v41,v43,v44),color=(0,0,1),material=frictMat)
  facets.extend((f11,f12,f21,f22,f31,f32,f41,f42))
O.bodies.append(facets)
mass = O.bodies[0].state.mass
for f in facets:
 f.state.mass = mass
 f.shape.color=(200,255,240)
 #f.state.blockedDOFs = 'XYZz'

# apply prestress to facets
#def addForces():
preStress = -3e9
for f in facets:
 n = f.shape.normal
 a = f.shape.area
 O.forces.addF(f.id,preStress*a*n)
 #O.forces.setPermF(f.id,preStress*a*n)
 #
print "force"

O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  #Bo1_Box_Aabb(),
  Bo1_Sphere_Aabb(),
  Bo1_GridConnection_Aabb(),
  Bo1_Facet_Aabb(),
 ]),
 InteractionLoop([
  Ig2_Sphere_Sphere_ScGeom(),
  #Ig2_Box_Sphere_ScGeom(),
  Ig2_GridNode_GridNode_GridNodeGeom6D(),
  Ig2_Sphere_GridConnection_ScGridCoGeom(),
  Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
  Ig2_Facet_Sphere_ScGeom(),
 ],
 [
  Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False), # internal cylinder physics
  Ip2_FrictMat_FrictMat_FrictPhys() # physics for external interactions, i.e., cylinder-cylinder, sphere-sphere, cylinder-sphere
  #Ip2_FrictMat_CpmMat_FrictPhys(),
 ],
 [
  Law2_ScGeom_FrictPhys_CundallStrack(), # contact law for sphere-sphere
  Law2_ScGridCoGeom_FrictPhys_CundallStrack(), # contact law for cylinder-sphere
  Law2_ScGeom6D_CohFrictPhys_CohesionMoment(), # contact law for "internal" cylinder forces
  Law2_GridCoGridCoGeom_FrictPhys_CundallStrack() # contact law for cylinder-cylinder interaction
  #Law2_ScGeom_CpmPhys_Cpm(),
 ]
 ),
 GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.1,label='ts'),
 #triax,
 #TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.1,label='newton'),
 #PyRunner(iterPeriod=1,command="addForces()"),
 PyRunner(command='plotAddData()',iterPeriod=10),
 #PyRunner(iterPeriod=20,command='history()',label='recorder'),
]

facet1 = facets[0] # or any other valid index
print(facet1.state.pos)
O.step()
print(facet1.state.pos) # should be different from the first one
#####################################################

Best Jan Stránský (honzik) said : #5

1) facets are fixed (state.blockedDOFs='xyzXYZ') by default.
You have to "unblock" some DOFs. In your last code, uncomment the line
f.state.blockedDOFs = 'XYZz'
and the position is different. Specifically it is (NaN,NaN,0) :-)

2) Further, you have to set non-zero mass to get non-NaN position

3)
addF add force only for one time step (="only a little")
setPermF set permanent force (you can change its value - also to 0 - during the simulation)

Your code is a bit complex for testing (e.g. interactions may play some role), use a real MWE for testing:
###
f = facet(((1,0,0),(0,1,0),(0,0,1)))
print("blockedDOFs",f.state.blockedDOFs)
f.state.blockedDOFs = "XYZ" # (!)
f.state.mass = 1 # (!)

O.bodies.append(f)
a,n = f.shape.area,f.shape.normal
p = -1
print("f",p*a*n)

print()
print("trying addF ... ")
O.forces.addF(0,p*a*n)
print("displ step",O.iter,f.state.displ())
O.step()
print("displ step",O.iter,f.state.displ())
O.step()
print("displ step",O.iter,f.state.displ())
O.step()
print("displ step",O.iter,f.state.displ())
print("addF applies 'constant velocity'")

print()
print("trying setPermF ... ")
O.forces.setPermF(0,p*a*n)
O.step()
print("displ step",O.iter,f.state.displ())
O.step()
print("displ step",O.iter,f.state.displ())
O.step()
print("displ step",O.iter,f.state.displ())
O.step()
print("displ step",O.iter,f.state.displ())
print("setPermF applies 'constant acceleration'")
###

cheers
Jan

onyourself (xxxe) said : #6

thanks, jan, that solved my problem.

onyourself (xxxe) said : #7

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