addForce didnt work for facet

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

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'),
]

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

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

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
#####################################################

thanks, jan, that solved my problem.

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