Commit 7d1cf1b1 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

hydro quebec 2d : lmgc

parent 87d5c1f6
from pylmgc90 import chipy
import numpy as np
class LmgcInterface(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(2)
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()
self._nbDisk = chipy.DISKx_GetNbDISKx()
self._d2bMap = chipy.DISKx_GetPtrDISKx2RBDY2()
self._position = np.empty([self._nbDisk,3], 'd')
self._velocity = np.empty([self._nbDisk,3], 'd')
self._externalF= np.empty([self._nbDisk,3], 'd')
self._volume = np.empty([self._nbDisk,1], 'd')
self._mass = np.empty([self._nbDisk,1], 'd')
for i in range(self._nbDisk):
self._volume[i] = np.pi * chipy.DISKx_GetContactorRadius(i+1)**2 * 2./3
self._mass[i] = chipy.RBDY2_GetBodyMass(int(self._d2bMap[i,0]))
def setBoundaryVelocity(self, tag, v) :
print("setBoundaryVelocity not implemented")
def write(self, odir, i, t) :
print("write not implemented")
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()
for i in range(self._nbDisk):
chipy.RBDY2_PutBodyVector('Fext_', int(self._d2bMap[i,0]), self._externalF[i,:])
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 the particles """
return self._volume
def mass(self):
""" Return the mass of the particles """
return self._mass
def position(self):
""" Get current position of contactors and return it """
for i in range(self._nbDisk):
self._position[i,:] = chipy.DISKx_GetContactorCoor(i+1)
self._position[i,2] = 0.
return self._position[:,:2]
def velocity(self):
""" Get current velocity of body of contactor and return it
Beware : it should not work very well with clusters !
"""
for i in range(self._nbDisk):
self._velocity[i,:] = chipy.RBDY2_GetBodyVector('V____',int(self._d2bMap[i,0]))
self._velocity[i,2] = 0.
return self._velocity[:,:2]
def externalForces(self):
""" Get an external forces array
"""
return self._externalF[:,:2]
def __del__(self):
"""
"""
chipy.ClosePostproFiles()
chipy.CloseDisplayFiles()
chipy.Finalize()
import os
import numpy as np
from pylmgc90 import pre_lmgc
import scontact2
datbox_path = 'DATBOX'
if not os.path.isdir(datbox_path):
os.mkdir(datbox_path)
sc = scontact2.ParticleProblem()
sc.load("deposited")
pillar_tag = 3
space_dim = 2
fric = 0
pb = []
x = sc.position()
v = sc.velocity()
v[:] = 0.
r = sc.particles()[:, 0]*0.9
pb = sc.segments()[sc.segmentTag() == 3,:][:,:2*space_dim].reshape(-1, 2)
rho_particle = np.mean(sc.particles()[:, 1]/(np.pi * r**2))
x_extremum = []
x_extremum.append(x[:,0].min()+r.min())
x_extremum.append(x[:,0].max()-r.min())
x_extremum.append(x[:,1].min()+r.min())
x_extremum = np.array(x_extremum)
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='Rxx2D',dimension=space_dim)
mat.addMaterial(mater)
mod.addModel(model)
for i in range(r.size) :
P = pre_lmgc.rigidDisk( r=r[i], center=x[i], model=model, material=mater, color='INxxx')
P.imposeInitValue(component=[1,2],value=list(v[i,:]))
avs.addAvatar(P)
r_min = np.min(r)
x_min = np.min(x[:,0])
x_max = np.max(x[:,0])
y_min = np.min(x[:,1])
y_max = np.max(x[:,1])
left = pre_lmgc.rigidJonc(axe1=y_max-y_min, axe2=r_min, center=[x_min-2*r_min,y_max],model=model,material=mater,color='OUTxx')
right= pre_lmgc.rigidJonc(axe1=y_max-y_min, axe2=r_min, center=[x_max+2*r_min,y_max],model=model,material=mater,color='OUTxx')
down = pre_lmgc.rigidJonc(axe1=(x_max-x_min)/2., axe2=r_min, center=[(x_max+x_min)/2., y_min-2*r_min],model=model,material=mater,color='OUTxx')
left.rotate(psi=-np.pi/2., center=left.nodes[1].coor)
right.rotate(psi=np.pi/2., center=right.nodes[1].coor)
left.imposeDrivenDof(component=[1,2,3],dofty='vlocy')
right.imposeDrivenDof(component=[1,2,3],dofty='vlocy')
down.imposeDrivenDof(component=[1,2,3],dofty='vlocy')
avs.addAvatar(left)
avs.addAvatar(right)
avs.addAvatar(down)
xmin = np.min(pb[:, 0])
xmax = np.max(pb[:, 0])
ymin = np.min(pb[:, 1])
ymax = np.max(pb[:, 1])
pillar = pre_lmgc.rigidPolygon(model,mater,np.zeros([space_dim]),color='PILAR',generation_type="full",vertices=np.array([(xmin,ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax)]))
avs.addAvatar(pillar)
#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,3],dofty='vlocy',ct=0)
pillar.imposeDrivenDof(component=[2],dofty='vlocy',ct=0)
#pillar.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 ="RBDY2", candidat ="DISKx", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY2", antagoniste="DISKx", colorAntagoniste='INxxx',
behav=clb_fric, alert=r.min())
OUT = pre_lmgc.see_table(CorpsCandidat ="RBDY2", candidat ="DISKx", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY2", antagoniste="JONCx", colorAntagoniste='OUTxx',
behav=clb_fric, alert=r.min())
PIL = pre_lmgc.see_table(CorpsCandidat ="RBDY2", candidat ="DISKx", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY2", antagoniste="POLYG", colorAntagoniste='PILAR',
behav=clb_fric, alert=r.min())
svs += IN
svs += OUT
svs += PIL
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)
pre_lmgc.writeTactBehav(tac,svs,chemin=datbox_path)
pre_lmgc.writeVlocRlocIni(chemin=datbox_path)
pre_lmgc.writePostpro(commands=post, parts=avs, path=datbox_path)
pre_lmgc.visuAvatars(avs)
import numpy as np
import shutil
class ScontactInterface :
def __init__(self, inputfile) :
import scontact2
p = scontact2.ParticleProblem()
self._p = p
p.load(inputfile)
p.velocity()[:, :] = 0;
self._r = np.mean(p.particles()[:, 0])
self._volume = np.pi * p.particles()[:, [0]]**2 * 2./3 ## 2./3 = 2d-3d correction
self._mass = p.particles()[:, [1]]
self._position = p.position()
self._velocity = p.velocity()
def volume(self):
return self._volume
def mass(self):
return self._mass
def position(self):
return self._position
def velocity(self):
return self._velocity
def setBoundaryVelocity(self, tag, v) :
p = self._p
p.disks()[p.diskTag() == tag, 2:4] = v
p.segments()[p.segmentTag() == tag, 4:6] = v
def iterate(self, dt, forces) :
self._p.externalForces()[:] = forces
self._p.iterate(self._r, dt, 1e-5)
def write(self, odir, i, t) :
outputname = "%s/part-%05i" % (odir, i)
self._p.write(outputname + "_tmp")
shutil.move(outputname + "_tmp", outputname)
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