# How to assign material properties？

Hi,
I have two scripts. They are used to generate a same sample and the difference of them is the way to assign properties for the material (sphere). The first program work properly but the second one dosen't.
I want to know what’s the matter with the second script.
Thanks.

The scripts are given below.

## 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 William (qfxx-123) said on 2022-05-13: #1

the first script

#########################################
### Defining parameters and variables ###
#########################################

### Matrix
Density1 = 2600
PoissonRatio1 = 0.5
frictionangle1 = 10
Young1 = 10e7
Num1 = 9000
Damp = 0.25

#Wall constants
WDensity = 0
WFrictionAngle = 0.0
WPoissonRatio = 0.5
WYoung = 50e9

#Packing variables
Size = 0.05
mn = Vector3(0, 0, 0)
mx = Vector3(Size, Size, Size)

#Confining variables
ConfPress = -190e3
ConfPress1 = ConfPress*0.9 # Pre-Compressure

########################################
#import necessary packages
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab

########################################
### Sample Preparing ##################################################################################################################################################################################
########################################

#Create materials for spheres and plates
WallMat = O.materials.append(FrictMat(young = WYoung , poisson = WPoissonRatio , frictionAngle = radians(WFrictionAngle) , density = WDensity))
wallIds = O.bodies.append(aabbWalls([mn,mx] , thickness = 0.001 , material = WallMat))

sp = pack.SpherePack()
sp.makeCloud(minCorner=mn, maxCorner=mx, rMean=0.0001, rRelFuzz=0, num=Num1)
sp.toSimulation(material = SphereMat)

O.usesTimeStepper=True

O.trackEnerty=True

#################################
##### Defining triaxil engines ##
#################################

###First step: Compression#######
triax1 = TriaxialStressController(
thickness = 0.001,
internalCompaction = False,
computeStressStrainInterval = 10,
goal1 = ConfPress1,
goal2 = ConfPress1,
goal3 = ConfPress1,
)

newton = NewtonIntegrator(damping = Damp)

### engine ###
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
GlobalStiffnessTimeStepper(active = 1 , timeStepUpdateInterval = 100 , timestepSafetyCoefficient = 0.8 , defaultDt = 4*utils.PWaveTimeStep()),
triax1,
newton,
PyRunner(realPeriod = 10 , command = 'checkUnbalanced()' , label = 'check'),
]

# Simulation stop conditions defination
def checkUnbalanced():
unb = unbalancedForce()
mStress = (triax1.stress(triax1.wall_right_id) + triax1.stress(triax1.wall_top_id) + triax1.stress(triax1.wall_front_id))/3

if unb<0.001 and abs(ConfPress1-mStress)/abs(ConfPress1)<0.01:
print('----------1 Finished----------')
print('----------Time:'+ str((time.time()-time0)/60) + 'min----------')
O.pause()

 Revision history for this message William (qfxx-123) said on 2022-05-13: #2

the second script

#########################################
### Defining parameters and variables ####
#########################################

### Matrix material
Density1 = 2600
PoissonRatio1 = 0.5
frictionangle1 = 10
Young1 = 10e7
Num1 = 9000
Damp = 0.25

#Wall constants
WDensity = 0
WFrictionAngle = 0.0
WPoissonRatio = 0.5
WYoung = 50e9

#Packing variables
Size = 0.05
mn = Vector3(0, 0, 0)
mx = Vector3(Size, Size, Size)

#Confining variables
ConfPress = -190e3
ConfPress1 = ConfPress*0.9 # Pre-Compressure

########################################
#import necessary packages
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab

########################################
### Sample Preparing ###################
########################################

#Create materials for spheres and plates
WallMat = O.materials.append(FrictMat(young = WYoung , poisson = WPoissonRatio , frictionAngle = radians(WFrictionAngle) , density = WDensity))
wallIds = O.bodies.append(aabbWalls([mn,mx] , thickness = 0.001 , material = WallMat))

sp = pack.SpherePack()
sp.makeCloud(minCorner=mn, maxCorner=mx, rMean=0.0001, rRelFuzz=0, num=Num1)
sp.toSimulation()

for b in O.bodies:
if isinstance(b.shape,Sphere):
b.shape.color = (0, 0, 1)
b.mat=O.materials[mat1]

O.usesTimeStepper=True

O.trackEnerty=True

#################################
##### Defining triaxil engines ##
#################################

###First step: Compression#######
triax1 = TriaxialStressController(
thickness = 0.001,
internalCompaction = False,
computeStressStrainInterval = 10,
goal1 = ConfPress1,
goal2 = ConfPress1,
goal3 = ConfPress1,
)

newton = NewtonIntegrator(damping = Damp)

### engine ###
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
GlobalStiffnessTimeStepper(active = 1 , timeStepUpdateInterval = 100 , timestepSafetyCoefficient = 0.8 , defaultDt = 4*utils.PWaveTimeStep()),
triax1,
newton,
PyRunner(realPeriod = 10 , command = 'checkUnbalanced()' , label = 'check'),
]

# Simulation stop conditions defination
def checkUnbalanced():
unb = unbalancedForce()
mStress = (triax1.stress(triax1.wall_right_id) + triax1.stress(triax1.wall_top_id) + triax1.stress(triax1.wall_front_id))/3
if unb<0.001 and abs(ConfPress1-mStress)/abs(ConfPress1)<0.01:
print('----------1 Finished----------')
print('----------Time:'+ str((time.time()-time0)/60) + 'min----------')
O.pause()

 Revision history for this message Jan Stránský (honzik) said on 2022-05-13: #3

Hello,

> The first program work properly but the second one dosen't.

please be more specific, what are the symptoms of "doesn't work properly".
Some errors? What errors?
Behavior is different? How?
... ?

Cheers
Jan

 Revision history for this message William (qfxx-123) said on 2022-05-13: #4

Sorry for my vague expression.

The 3D model in primary view will disappear suddenly when the second script is run.

 Revision history for this message  Jan Stránský (honzik) said on 2022-05-13: #5

Thanks, now it is clear.

In the second script, you call sp.toSimulation() without material argument.
Therefore the last O.materials.appended material is used, i.e. WallMat.
WallMat has zero density. Therefore spheres has zero mass.
Therefore they disappear, because from Newton's law a=F/m zero mass leads to infinite acceleration and then to NaN position and are not displayed.

The code
for b in O.bodies:
....
b.mat=O.materials[mat1]

has no effect on already computed values (like mass), you would have to update them "manually".

Cheers
Jan

 Revision history for this message William (qfxx-123) said on 2022-05-14: #6

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