# addForce didnt work for facet

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.state.mass
for f in facets:
f.state.mass = mass
f.shape.color=(200,255,240)
f.state.blockedDOFs = 'XYZz'

for f in facets:
n = f.shape.normal
a = f.shape.area
##########################
my whole code is a bit long so i upload related parts of it, thanks.

## Question information

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

Hello,

1) please provide a MWE , 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:
- 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 , 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

 onyourself (xxxe) said on 2020-10-15: #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 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
# 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,
)

#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.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):
x=b.state.pos
y=b.state.pos
z=b.state.pos
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,c)>=min(x0,x1)
and max(x0,x1)>=min(c,c)
and max(c,c)>=min(z0,z1)
and max(z0,z1)>=min(c,c)):
x12=(x1-x0)
z12=(z1-z0)
x13=(c-x0)
z13=(c-z0)
x14=(c-x0)
z14=(c-z0)
cross1=x12*z13-x13*z12
cross2=x12*z14-x14*z12
x34=c-c
z34=c-c
#x23=c-x1
#z23=c-z1
x32=x1-c
z32=z1-c
#x24=c-x1
#z24=c-z1
x31=x0-c
z31=z0-c
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):
if r>0.012e-3:
x=b.state.pos
y=b.state.pos
z=b.state.pos
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.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
for f in facets:
n = f.shape.normal
a = f.shape.area
'''
for f in facets:
n = f.shape.normal
a = f.shape.area
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 on 2020-10-15: #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
###

f = facets # 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 on 2020-10-16: #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 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.state.mass
for f in facets:
f.state.mass = mass
f.shape.color=(200,255,240)
#f.state.blockedDOFs = 'XYZz'

# apply prestress to facets
preStress = -3e9
for f in facets:
n = f.shape.normal
a = f.shape.area
#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=20,command='history()',label='recorder'),
]

facet1 = facets # or any other valid index
print(facet1.state.pos)
O.step()
print(facet1.state.pos) # should be different from the first one
##################################################### Jan Stránský (honzik) said on 2020-10-16: #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)
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("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()
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 on 2020-10-21: #6

thanks, jan, that solved my problem.

 onyourself (xxxe) said on 2020-10-21: #7

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