Commit a2169f9e authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
parents e01591ec abada195
import numpy as np
from pylmgc90 import chipy
class LmgcParticles(object):
"""
A class definition holding data of lmgc90
in a compatible form to use coupling with
gmsh
Attributs are :
- number of disks
- volume of diskx
- reference on diskx2rbdy2 map
The methods are :
- iterate() : do one time iteration
- position() : give back position of all diskx
- velocity() : give back position of all diskx
- externalForces() : update external forces on particles
- getMeanRadius() : return the mean radius of particles
- getVolume() : give back the volume of all diskx
- writeHeader() : write header file for plot
"""
def __init__(self,space_dim, dt, theta) :
"""
Initialize LMGC90
"""
### Set dimension in chipy for dummies ###
chipy.SetDimension(space_dim)
chipy.Initialize()
### definition des parametres du calcul ###
chipy.utilities_logMes('INIT TIME STEPPING')
chipy.TimeEvolution_SetTimeStep(dt)
chipy.Integrator_InitTheta(theta)
### lecture du modele ###
chipy.utilities_logMes('READ BEHAVIOURS')
chipy.ReadBehaviours()
chipy.utilities_logMes('READ BODIES')
chipy.ReadBodies()
chipy.utilities_logMes('LOAD BEHAVIOURS')
chipy.LoadBehaviours()
chipy.utilities_logMes('READ INI DOF')
chipy.ReadIniDof()
chipy.utilities_logMes('READ DRIVEN DOF')
chipy.ReadDrivenDof()
chipy.utilities_logMes('LOAD TACTORS')
chipy.LoadTactors()
chipy.utilities_logMes('READ INI Vloc Rloc')
chipy.ReadIniVlocRloc()
##### ecriture paranoiaque du modele ###
##utilities_logMes('WRITE BODIES')
##WriteBodies()
##utilities_logMes('WRITE BEHAVIOURS')
##WriteBehaviours()
##utilities_logMes('WRITE DRIVEN DOF')
##WriteDrivenDof()
chipy.OpenPostproFiles()
chipy.OpenDisplayFiles()
chipy.WriteDisplayFiles()
chipy.ComputeMass()
# list ignored spheres in gmsh
self._s2s = []
self._space_dim = space_dim
if space_dim == 2:
self._nbDisk = chipy.DISKx_GetNbDISKx()
self._d2bMap = chipy.DISKx_GetPtrDISKx2RBDY2()
self._getCoor = chipy.DISKx_GetContactorCoor
self._getVect = chipy.RBDY2_GetBodyVector
self._setVect = chipy.RBDY2_PutBodyVector
nb_dof = 3
for i in xrange(shper.nb_Disk):
if chipy.RBDY2_GetContactorColor(int(self_d2bMap[i,0]),int(self._d2bMap[i,1])) == "PILAR":
self._s2s.append(i)
else :
self._nbDisk = chipy.SPHER_GetNbSPHER()
self._d2bMap = chipy.SPHER_GetPtrSPHER2RBDY3()
self._getCoor = chipy.SPHER_GetContactorCoor
self._getVect = chipy.RBDY3_GetBodyVector
self._setVect = chipy.RBDY3_PutBodyVector
nb_dof = 6
for i in xrange(self._nbDisk):
if chipy.RBDY3_GetContactorColor(int(self._d2bMap[i,0]),int(self._d2bMap[i,1])) == "PILAR":
self._s2s.append(i)
nb_skip = len(self._s2s)
self._position = np.empty([self._nbDisk-nb_skip,nb_dof], 'd')
self._velocity = np.empty([self._nbDisk-nb_skip,nb_dof], 'd')
self._externalF= np.empty([self._nbDisk-nb_skip,nb_dof], 'd')
self._volume = np.empty([self._nbDisk-nb_skip,1], 'd')
if space_dim == 2:
ishift = 0
for i in xrange(self._nbDisk):
if i in self._s2s:
ishift += 1
continue
self._volume[i-ishift] = np.pi * chipy.DISKx_GetContactorRadius(i+1)**2
else:
ishift = 0
for i in xrange(self._nbDisk):
if i in self._s2s:
ishift += 1
continue
self._volume[i-ishift] = 4. * np.pi * chipy.SPHER_GetContactorRadius(i+1)**3 / 3.
def iterate(self,freq_detect,freq_write,freq_display,ref_radius,**solver_args):
"""
Do one step of a LMGC90 computation.
Last argument is a dictionnary holding nlgs_solver list of arguments.
"""
#
chipy.IncrementStep()
chipy.ComputeFext()
ishift = 0
for i in xrange(self._nbDisk):
if i in self._s2s:
ishift += 1
continue
self._setVect('Fext_', int(self._d2bMap[i,0]), self._externalF[i-ishift,:])
chipy.ComputeBulk()
chipy.ComputeFreeVelocity()
chipy.SelectProxTactors(freq_detect)
chipy.RecupRloc()
chipy.ExSolver(**solver_args)
chipy.StockRloc()
chipy.ComputeDof()
if freq_write > 0:
chipy.WriteOutDof(freq_write)
chipy.WriteOutVlocRloc(freq_write)
else:
chipy.WriteLastDof()
chipy.WriteLastVlocRloc()
chipy.WritePostproFiles()
if freq_display > 0:
chipy.WriteDisplayFiles(freq_display,ref_radius)
chipy.UpdateStep()
def volume(self):
""" Return the volume of particles """
return self._volume
def position(self):
""" Get current position of contactors and return it """
ishift = 0
for i in xrange(self._nbDisk):
if i in self._s2s:
ishift += 1
continue
self._position[i-ishift,:] = self._getCoor(i+1)
self._position[i-ishift,self._space_dim] = 0.
return self._position[:,:self._space_dim]
def velocity(self):
""" Get current velocity of body of contactor and return it
Beware : it should not work very well with clusters !
"""
ishift = 0
for i in xrange(self._nbDisk):
if i in self._s2s:
ishift += 1
continue
self._velocity[i-ishift,:] = self._getVect('V____',int(self._d2bMap[i,0]))
self._velocity[i-ishift,self._space_dim:] = 0.
return self._velocity[:,:self._space_dim]
def externalForces(self):
""" Get an external forces array
"""
return self._externalF[:,:self._space_dim]
def __del__(self):
"""
"""
chipy.ClosePostproFiles()
chipy.CloseDisplayFiles()
chipy.Finalize()
import numpy as np
import os, time
from dgpy import *
import runf
outputdir ="loutput_run"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
class meshVelocity :
def __init__(self) :
self.Y = 0
self.V = 0
def cb(self, cm, v, x) :
l = 0.1
L = 1
h = 0.5
H = 1.5
R = np.hypot(x[:, 0], x[:, 1])
fx = np.minimum(1 - (R - l) / (L - l), 1.)
HX = self.Y * fx
HH = np.where(x[:, 2] > HX, H - HX , H + HX)
fy = np.minimum(1 - (np.abs(x[:, 2] - HX) - h) / (HH - h), 1.)
v[:, 0:2] = 0
v[:, 2] = fx * fy * self.V
def update(self, dt) :
self.Y += self.V * dt
tEnd = 3.
tCycle = 3.
fric = 0.
rho_particle = 2.
g = 1.
gravity = [0., 0., -g]
R_cyl = 1.
from pylmgc90 import pre_lmgc
datbox_path = 'DATBOX'
if not os.path.isdir(datbox_path):
os.mkdir(datbox_path)
f = open('deposited')
pillar_tag = 3
l = f.readline()
space_dim = int(l.split()[1])
x = []
v = []
r = []
pb = []
pb_n = []
pb_s = []
ipb = 0
for l in f.readlines():
if l[0] == 'P':
l = l.split()
x.append( map(float,l[1:space_dim+1]) )
v.append( map(float,l[space_dim+1:2*space_dim+1]) )
r.append( float(l[2*space_dim+1]) )
elif l[0] == 'T':
l = l.split()
if int(l[-1]) == pillar_tag:
s = []
for i in xrange(3):
pTi = map(float,l[i*space_dim+1:(i+1)*space_dim+1])
if pTi not in pb_n:
pb_n.append(pTi)
ipb += 1
pb.append(ipb)
s.append(ipb)
else:
for k, n in enumerate(pb_n):
if n==pTi:
s.append(k+1)
break
pb_s.append(s)
x = np.array(x)
r = np.array(r)
pb_n = np.array(pb_n)
pb_s = np.array(pb_s)
avs = pre_lmgc.avatars()
mat = pre_lmgc.materials()
mod = pre_lmgc.models()
svs = pre_lmgc.see_tables()
tac = pre_lmgc.tact_behavs()
mater = pre_lmgc.material(name='STONE',type='RIGID',density=rho_particle)
model = pre_lmgc.model(name='rigid', type='MECAx', element='Rxx3D',dimension=space_dim)
mat.addMaterial(mater)
mod.addModel(model)
for i in xrange(r.size) :
P = pre_lmgc.rigidSphere( r=r[i], center=x[i], model=model, material=mater, color='INxxx')
#P.imposeInitValue(component=[1,2,3],value=v[i])
avs.addAvatar(P)
r_min = np.min(r)
r_max = np.max(r)
z_max = np.max(x[:,2]+r)
z_min = np.min(x[:,2]-r)
lyc = pre_lmgc.rigidCylinder( R_cyl, 1.5*(z_max-z_min), [0.,0.,(3.*z_max+z_min)/4.], model, mater, 'OUTxx', is_Hollow=True )
down = pre_lmgc.rigidPlan(axe1=R_cyl, axe2=R_cyl, axe3=r_min, center=[0., 0., z_min-r_min], model=model, material=mater, color='OUTxx')
lyc.imposeDrivenDof(component=range(1,7),dofty='vlocy')
down.imposeDrivenDof(component=range(1,7),dofty='vlocy')
avs.addAvatar(lyc)
avs.addAvatar(down)
pillar = pre_lmgc.rigidPolyhedron(model,mater,np.zeros([space_dim]),color='PILAR',generation_type="full",vertices=pb_n,faces=pb_s)
#avs.addAvatar(pillar)
h_pil = np.max(pb_n[:,2]) - np.min(pb_n[:,2])
z_pil = ( np.max(pb_n[:,2]) + np.min(pb_n[:,2]) ) / 2.
r_pil = np.max( np.linalg.norm( pb_n[:,:2], axis=1 ) )
pil = pre_lmgc.rigidCylinder( r_pil, h_pil, [0.,0.,z_pil], model, mater, 'PILAR' )
pil.addContactors( 'SPHER','PILAR','all', byrd=r_pil, shift=[0.,0.,-h_pil/2.] )
pil.computeRigidProperties()
avs.addAvatar(pil)
#sys.exit()
def imposedVelocity(t):
tCycle = 3.
return (-1 if ((t % tCycle) < (tCycle / 2)) else 1) * 1/ (tCycle/2)
times = [0.,tCycle/2.-1.e-7,tCycle/2.]
for i in xrange(int(2*tEnd/tCycle)):
times.extend(map(lambda x: x+tCycle/2.,times[-2:]))
pre_lmgc.writeEvolution(f=imposedVelocity, instants=times, path='DATBOX/', name='v_pillar.dat')
pillar.imposeDrivenDof(component=[1,2,4,5,6],dofty='vlocy')
pillar.imposeDrivenDof(type='evolution',component=space_dim,dofty='vlocy',evolutionFile='v_pillar.dat')
pil.imposeDrivenDof(component=[1,2,4,5,6],dofty='vlocy')
pil.imposeDrivenDof(type='evolution',component=space_dim,dofty='vlocy',evolutionFile='v_pillar.dat')
clb_fric = pre_lmgc.tact_behav(name='iqsc0',type='IQS_CLB',fric=fric)
tac += clb_fric
IN = pre_lmgc.see_table(CorpsCandidat ="RBDY3", candidat ="SPHER", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY3", antagoniste="SPHER", colorAntagoniste='INxxx',
behav=clb_fric, alert=r.max())
OUT = pre_lmgc.see_table(CorpsCandidat ="RBDY3", candidat ="SPHER", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY3", antagoniste="PLANx", colorAntagoniste='OUTxx',
behav=clb_fric, alert=r.min())
CYL = pre_lmgc.see_table(CorpsCandidat ="RBDY3", candidat ="SPHER", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY3", antagoniste="DNLYC", colorAntagoniste='OUTxx',
behav=clb_fric, alert=r.max())
PIL = pre_lmgc.see_table(CorpsCandidat ="RBDY3", candidat ="SPHER", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY3", antagoniste="POLYR", colorAntagoniste='PILAR',
behav=clb_fric, alert=2.*r.max())
PIL2= pre_lmgc.see_table(CorpsCandidat ="RBDY3", candidat ="SPHER", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY3", antagoniste="CYLND", colorAntagoniste='PILAR',
behav=clb_fric, alert=2.*r.max())
PIL3= pre_lmgc.see_table(CorpsCandidat ="RBDY3", candidat ="SPHER", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY3", antagoniste="SPHER", colorAntagoniste='PILAR',
behav=clb_fric, alert=2.*r.max())
svs += IN
svs += CYL
svs += OUT
svs += PIL
svs += PIL2
svs += PIL3
post = pre_lmgc.postpro_commands()
solv = pre_lmgc.postpro_command(type='SOLVER INFORMATIONS', step=1)
post.addCommand(solv)
# file writting
pre_lmgc.writeBodies(avs,chemin=datbox_path)
pre_lmgc.writeDrvDof(avs,chemin=datbox_path)
pre_lmgc.writeDofIni(avs,chemin=datbox_path)
pre_lmgc.writeModels(mod,chemin=datbox_path)
pre_lmgc.writeBulkBehav(mat,chemin=datbox_path,gravy=gravity)
pre_lmgc.writeTactBehav(tac,svs,chemin=datbox_path)
pre_lmgc.writeVlocRlocIni(chemin=datbox_path)
pre_lmgc.writePostpro(commands=post, parts=avs, path=datbox_path)
from lmgc_particles import *
t = 0.
dt = 0.005 * tCycle
theta = 0.5
freq_write = 5
freq_display = 5
ref_radius = 0.1e-1
freq_detect = 1
solver_params = { 'type' :'Stored_Delassus_Loops ',
'norm' : 'Quad ',
'conv' : 1e-5,
'relax' : 1.0,
'gsit1' : 20,
'gsit2' : 2000
}
#chipy.nlgs_SetWithQuickScramble()
chipy.checkDirectories()
particles = LmgcParticles(space_dim,dt,theta)
outf = freq_display
mv = meshVelocity()
hydro = runf.hydro("mesh.msh", mv.cb, dt, np.mean(r), rhop=rho_particle, g=gravity[:space_dim])
mv.V = 0.5/dt
hydro.moveMesh()
mv.update(dt)
t = 0
ii = 0
outputname = "%s/part-%05i" % (outputdir, 0)
tic = time.clock()
while t < tEnd :
mv.V = (-1 if ((t % tCycle) < (tCycle / 2)) else 1) * 1/ (tCycle/2)
hydro.moveMesh()
print("solve")
hydro.solve(particles.volume(), particles.position(), particles.velocity(), particles.externalForces())
print("solved")
particles.iterate(freq_detect,freq_write,freq_display,ref_radius,**solver_params)
mv.update(dt)
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf)
hydro.writeSolution(outputdir, ioutput, t, filetype="vtk")
ii += 1
print("%.2g/%.2g" % (t, tEnd))
print("cpu time : %g" % (time.clock() - tic))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment