Question #699471 on Yade changed:
https://answers.launchpad.net/yade/+question/699471
Mithushan Soundaranathan posted a new comment:
Here is code for saved simulation, which loaded into the new simulation:
from __future__ import print_function
from yade import utils, plot, timing
from yade import pack
import pandas as pd
import numpy as np
from PIL import Image
from yade import pack, export
from scipy.interpolate import interp1d
readParamsFromTable(k1=2e5,kp=2e6)
from yade.params.table import *
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import os
import csv
from matplotlib.pyplot import figure
from scipy.interpolate import interp1d
from pylab import *
from scipy.optimize import curve_fit
o = Omega()
save=1
# Physical parameters
fr = 0.41
rho = 1561
Diameter = 0.0012/2
D=Diameter
r1 = Diameter/2
#r2 = Diameter/2
k1 = 1.26e5
kp = 12*k1
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)
if save==1:
o.dt = 1.0e-7
else:
o.dt=1.0e-5
particleMass = 4.0/3.0*math.pi*r1*r1*r1*rho
Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
m_tot=0.0005
#no_p=m_tot/particleMass
no_p=1000
Tab_rad=0.005
Cyl_height=0.015
cross_area=math.pi*(Tab_rad**2)
Comp_press= 1.2e8
Comp_force=Comp_press*cross_area
##Single particle swelling model
def model(r,t,Q_max,rho_t,rho_w,r_0,D):
Q=((rho_w*(r**3))/(rho_t*(r_0**3)))-(rho_w/rho_t)+1;
drdt =((D*rho_t)/(r*rho_w))*((Q_max-Q)/Q);
return drdt
P=[1.45,rho,1000,396.39e-12]
#*************************************
# Add material
mat1 = o.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1,
kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))
# Spheres for compression and walls
sp=pack.SpherePack()
sp.makeCloud((-5.0*Diameter,-5.0*Diameter,-12*Diameter),(5.0*Diameter,5.0*Diameter,7.0*Diameter),
rMean=Diameter/2.0,rRelFuzz=0.18,num=1000)
sp.toSimulation(material=mat1)
walls=o.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=Cyl_height,segmentsNumber=20,wallMask=6,material=mat1))
r_save=[]
radius=[]
initial_save=[]
size_save=[]
radius.append(0)
swell_t=np.zeros((1,no_p))
i=0
for b in O.bodies:
if isinstance(b.shape, Sphere):
radius.append(b.shape.radius)
r_save.append(radius)
# Add engines
o.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
Bo1_Wall_Aabb(),
Bo1_Facet_Aabb()
]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
Ig2_Facet_Sphere_ScGeom(),
Ig2_Wall_Sphere_ScGeom()],
[Ip2_LudingMat_LudingMat_LudingPhys()],
[Law2_ScGeom_LudingPhys_Basic()]
),
NewtonIntegrator(damping=0.1, gravity=[0, 0, -9.81]),
PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
#DeformControl(label="DefControl")
]
def checkForce():
# at the very start, unbalanced force can be low as there is only few
# contacts, but it does not mean the packing is stable
if O.iter < 400000:
return
# the rest will be run only if unbalanced is < .1 (stabilized packing)
timing.reset()
if unbalancedForce() > 1:
return
# add plate at upper box side
highSphere = 0.0
for b in O.bodies:
if highSphere < b.state.pos[2] and isinstance(b.shape, Sphere):
highSphere = b.state.pos[2]
else:
pass
#o.bodies.append(wall(highSphere, axis=2, sense=-1, material=mat1))
#global plate
#plate = o.bodies[-1] # the last particles is the plat
#plate.state.vel = (0, 0, -2)
global upper_punch
upper_punch=O.bodies.append(geom.facetCylinder((0,0,((-Cyl_height/2)+0.0001)+utils.aabbDim()[2]),Tab_rad-.00001,0,segmentsNumber=50,wallMask=1))
for i in upper_punch:
body= O.bodies[i]
body.state.vel = (0,0,-2)
fCheck.command = 'unloadPlate()'
def unloadPlate():
force_up=0
for i in upper_punch:
body= O.bodies[i]
force_up=force_up+abs(O.forces.f(body.id)[2])
# if the force on plate exceeds maximum load, start unloading
# if abs(O.forces.f(plate.id)[2]) > 5e-2:
if force_up > Comp_force:
for i in upper_punch:
body= O.bodies[i]
body.state.vel = (0,0,2)
#plate.state.vel *= -1
# next time, do not call this function anymore, but the next one
# (stopUnloading) instead
fCheck.command = 'stopUnloading()'
def stopUnloading():
for i in upper_punch:
body=O.bodies[i]
pos_up=body.state.pos
if pos_up[2] > Cyl_height/2:
for i in upper_punch:
body= O.bodies[i]
body.state.vel = (0,0,2)
#plate.state.vel = (0, 0, 0)
initial_save.append(utils.aabbExtrema()[1][2]+r1)
initial_save.append(O.iter)
#o.engines = o.engines+[PyRunner(command='ParticleSwelling()',
iterPeriod=100000)]
#for b in O.bodies:
#if isinstance(b.shape,Wall): O.bodies.erase(b.id) # erase punch
for j in upper_punch: O.bodies.erase(j)
#o.pause()
#if save==1:
#utils.saveVars('somethingtest',plate=plate)
#o.save('Tablet_swelling.xml')
fCheck.command = 'ParticleSwelling()'
def ParticleSwelling():
if save==1:
utils.saveVars('somethingtest',plate=upper_punch)
o.save('Tablet_swelling.xml')
o.pause()
##time_current=(O.iter-initial_save[1])*o.dt
##Liq_pos=0.24*(time_current**0.574)
##Liq_pos=Liq_pos*1e-3 #convert to m
##radius=[]
##radius.append(time_current)
##for b in O.bodies:
##if isinstance(b.shape, Sphere):
##par_pos=(initial_save[0]-b.state.pos[2])
##k=b.id
##r_now=b.shape.radius
##if Liq_pos>=par_pos:
##r_0=r_save[0][k+1]
##if swell_t[0][k]==0:
##swell_t[0][k]=time_current
##radius.append(r_save[0][k+1])
##continue
##time=time_current-swell_t[0][k]
##t = np.linspace(0,time)
##r = odeint(model,r_0,t,args=(P[0],P[1],P[2],r_0,P[3],))
##r_new=float(r[-1])
##radius.append(r_new)
##b.shape.radius = float(r[-1])
##f=float(r[-1])/r_new
##mcurrent=b.state.mass
##mnew=mcurrent*(f*f*f)
##b.state.mass=mnew
##Icurrent=b.state.inertia
##Inew=Icurrent*(f*f*f*f*f)
##b.state.inertia=Inew
##elif Liq_pos<par_pos:
##radius.append(r_save[0][k+1])
##r_save.append(radius)
##size_current=[time_current,utils.aabbDim()[2]]
##size_save.append(size_current)
##if time_current>20:
##o.pause()
##radius_data=pd.DataFrame(r_save)
##path_save='/home/mithushan/Swelling'
##base_filename='PH101_swelling_data.csv'
##radius_data.to_csv(os.path.join(path_save,base_filename))
##size_data=pd.DataFrame(size_save, columns=['time','height'])
##base_filename='PH101_height.csv'
##size_data.to_csv(os.path.join(path_save,base_filename))
O.run();
O.wait()
#if save==1:
#o.save('Tablet_swelling.xml')
if save==0:
o.load('Tablet_swelling.xml')
utils.loadVars('somethingtest')
from yade.params.somethingtest import *
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
