Commit 9aee5105 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

hydro quebec 3d lmgc : almost there

parent f39056ba
......@@ -27,7 +27,7 @@ class LmgcInterface(object):
- writeHeader() : write header file for plot
"""
def __init__(self,space_dim, dt) :
def __init__(self) :
"""
Initialize LMGC90
"""
......@@ -192,7 +192,6 @@ def scontactToLmgc90(filename):
sc = scontact2.ParticleProblem()
sc.load(filename)
space_dim = 2
fric = 0
x = sc.position()
......@@ -208,7 +207,7 @@ def scontactToLmgc90(filename):
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)
model = pre_lmgc.model(name='rigid', type='MECAx', element='Rxx2D',dimension=2)
mat.addMaterial(mater)
mod.addModel(model)
......@@ -220,8 +219,8 @@ def scontactToLmgc90(filename):
seg = sc.segments()
bndtags = set()
for i in range(seg.shape[0]) :
x0 = np.array(seg[i, 0:space_dim])
x1 = np.array(seg[i, space_dim:2*space_dim])
x0 = np.array(seg[i, 0:2])
x1 = np.array(seg[i, 2:4])
tag = sc.getTagName(int(sc.segmentTag()[i]))
if len(tag) > 5 :
print("maximum tag name length in lmgc is 5, '%s' is too long" % tag)
......@@ -231,7 +230,7 @@ def scontactToLmgc90(filename):
t = (x0 - x1)/np.linalg.norm(x0 - x1)
n = np.array([-t[1], t[0]]) * 1e-12
vs = np.array([x0-n, x0 + n, x1 + n, x1 - n])
av = pre_lmgc.rigidPolygon(model,mater,np.zeros([space_dim]),color=tag,generation_type="full",vertices=vs)
av = pre_lmgc.rigidPolygon(model,mater,np.zeros([2]),color=tag,generation_type="full",vertices=vs)
av.imposeDrivenDof(component=[1,2,3],dofty='vlocy')
avs.addAvatar(av)
......
from pylmgc90 import chipy
import numpy as np
import os
class LmgcInterface(object):
"""
A class definition holding data of lmgc90
in a compatible form to use coupling with
gmsh
Attributs are :
- number of spheres
- volume of spheres
- reference on spher2rbdy3 map
The methods are :
- iterate() : do one time iteration
- position() : give back position of all spher
- velocity() : give back position of all spher
- externalForces() : update external forces on particles
- getMeanRadius() : return the mean radius of particles
- getVolume() : give back the volume of all spher
- writeHeader() : write header file for plot
"""
def __init__(self) :
"""
Initialize LMGC90
"""
chipy.checkDirectories()
### Set dimension in chipy for dummies ###
chipy.SetDimension(3)
chipy.Initialize()
### definition des parametres du calcul ###
### 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._nbSpher = chipy.SPHER_GetNbSPHER()
self._d2bMap = chipy.SPHER_GetPtrSPHER2RBDY3()
self._position = np.zeros([self._nbSpher,3], 'd')
self._velocity = np.zeros([self._nbSpher,6], 'd')
self._externalF= np.zeros([self._nbSpher,6], 'd')
self._volume = np.zeros([self._nbSpher,1], 'd')
self._mass = np.zeros([self._nbSpher,1], 'd')
for i in range(self._nbSpher):
self._volume[i] = np.pi * chipy.SPHER_GetContactorRadius(i+1)**3 * 4./3
self._mass[i] = chipy.RBDY3_GetBodyMass(int(self._d2bMap[i,0]))
self._tag2bnd = {}
self._tag2id = {}
for i in range(chipy.RBDY3_GetNbRBDY3()):
for t in range(chipy.RBDY3_GetNbContactor(i+1)):
tag = chipy.RBDY3_GetContactorColor(i+1,1).rstrip('_')
self._tag2id.setdefault(tag, []).append(i)
self._freq_detect = 1
self._solver_params = { 'type' :'Stored_Delassus_Loops ',
'norm' : 'Quad ',
'conv' : 1e-5,
'relax' : 1.,
'gsit1' : 200,
'gsit2' : 200
}
def setBoundaryVelocity(self, tag, v) :
self._tag2bnd[tag] = v
def write(self, odir, i, t) :
chipy.WriteOutDof(1)
chipy.WriteOutVlocRloc(1)
chipy.WritePostproFiles()
chipy.WriteDisplayFiles(1, 0.01)
#def iterate(self,freq_detect,freq_write,freq_display,ref_radius,**solver_args):
def iterate(self, dt, forces):
"""
Do one step of a LMGC90 computation.
"""
chipy.TimeEvolution_SetTimeStep(dt)
chipy.Integrator_InitTheta(1.)
self._externalF[:,:3] = forces
chipy.IncrementStep()
for tag, v in self._tag2bnd.iteritems() :
for i in self._tag2id[tag] :
for j, iv in enumerate(v):
chipy.RBDY3_SetVlocyDrivenDof(i+1,j+1,iv)
chipy.ComputeFext()
for i in range(self._nbSpher):
chipy.RBDY3_PutBodyVector('Fext_', int(self._d2bMap[i,0]), self._externalF[i,:])
chipy.ComputeBulk()
chipy.ComputeFreeVelocity()
chipy.SelectProxTactors(self._freq_detect)
chipy.RecupRloc()
chipy.ExSolver(**self._solver_params)
chipy.StockRloc()
chipy.ComputeDof()
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._nbSphere):
self._position[i,:] = chipy.SPHER_GetContactorCoor(i+1)
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._nbSPher):
self._velocity[i,:] = chipy.RBDY3_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()
def scontactToLmgc90(filename):
import scontact3
from pylmgc90 import pre_lmgc
datbox_path = 'DATBOX'
if not os.path.isdir(datbox_path):
os.mkdir(datbox_path)
sc = scontact3.ParticleProblem()
sc.load(filename)
fric = 0
x = sc.position()
v = sc.velocity()
v[:] = 0.
r = sc.particles()[:, 0]
rho_particle = np.mean(sc.particles()[:, 1]/(np.pi * r**3 * 4./3))
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=3)
mat.addMaterial(mater)
mod.addModel(model)
for i in range(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=list(v[i,:]))
avs.addAvatar(P)
tri = sc.triangles()
bndtags = set()
for i in range(tri.shape[0]) :
x0 = np.array(tri[i, 0:3])
x1 = np.array(tri[i, 3:6])
x2 = np.array(tri[i, 6:9])
tag = sc.getTagName(int(sc.triangleTag()[i]))
if len(tag) > 5 :
print("maximum tag name length in lmgc is 5, '%s' is too long" % tag)
exit()
tag += "_" * (5 - len(tag))
bndtags.add(tag)
t0 = (x0 - x1)/np.linalg.norm(x0 - x1)
t1 = (x0 - x2)/np.linalg.norm(x0 - x2)
n = np.cross(t0, t1)
n *= 1e-12 / np.linalg.norm(n)
vs = np.array([x0-n, x1 - n, x2 - n, x0 + n, x1 + n, x2 + n])
fs = np.array([[0, 2, 1], [3, 4, 5], [0, 1, 3], [3, 1, 4], [1, 5, 4], [1, 2, 5], [0, 3, 2], [3, 5, 2]]) + 1
av = pre_lmgc.rigidPolyhedron(model,mater,np.zeros([3]),color=tag,generation_type="full",vertices=vs, faces=fs)
av.imposeDrivenDof(component=[1,2,3],dofty='vlocy')
avs.addAvatar(av)
clb_fric = pre_lmgc.tact_behav(name='iqsc0',type='IQS_CLB',fric=fric)
tac += clb_fric
svs += pre_lmgc.see_table(CorpsCandidat ="RBDY3", candidat ="SPHER", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY3", antagoniste="SPHER", colorAntagoniste='INxxx',
behav=clb_fric, alert=r.min())
for tag in bndtags :
svs += pre_lmgc.see_table(CorpsCandidat ="RBDY3", candidat ="SPHER", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY3", antagoniste="POLYR", colorAntagoniste=tag,
behav=clb_fric, alert=r.min())
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)
......@@ -31,9 +31,9 @@ if not os.path.isdir(outputdir) :
mv = MeshVelocity()
# LMGC
from pyfefp import lmgcInterface
lmgcInterface.scontactToLmgc90("deposited")
p = lmgcInterface.LmgcInterface(2, dt)
from pyfefp import lmgc2Interface
lmgc2Interface.scontactToLmgc90("deposited")
p = lmgc2Interface.LmgcInterface()
# SCONTACT
#from pyfefp import scontactInterface
......
......@@ -33,13 +33,13 @@ if not os.path.isdir(outputdir) :
mv = MeshVelocity()
# LMGC
#from pyfefp import lmgcInterface
#lmgcInterface.scontactToLmgc90("deposited")
#p = lmgcInterface.LmgcInterface(3, dt)
from pyfefp import lmgc3Interface
lmgc3Interface.scontactToLmgc90("deposited")
p = lmgc3Interface.LmgcInterface()
# SCONTACT
from pyfefp import scontact3Interface
p = scontact3Interface.ScontactInterface("deposited")
#from pyfefp import scontact3Interface
#p = scontact3Interface.ScontactInterface("deposited")
rhop = np.mean(p.mass() / p.volume())
fluid = pyfefp.FluidSolver("mesh.msh", dt, rhop=rhop, rho=1, mu=1, g = [0, 0, -1], imex=False)
......
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