Commit 8e5cb2cd authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

merge 3d and 2d code

parent bb81513a
ALL: libmbfluid.so libscontact2.so libscontact3.so
CFLAGS=-Wno-unused-function -g3 -march=native -mtune=native -fPIC -std=gnu99
LDFLAGS=-shared -lm
FLUID_C=src/fluid_problem.c src/mesh.c src/mesh_find.c src/quadtree.c src/hxt_linear_system.c src/hxt_linear_system_lu.c src/hxt_message.c
FLUID_H=src/tools.h src/fluid_problem.h src/hxt_linear_system.h src/mesh_find.h src/mesh.h src/quadtree.h src/vector.h
USE_PETSC=1
ifeq ($(USE_PETSC),1)
FLUID_C += src/hxt_linear_system_petsc.c
CFLAGS += -I${PETSC_DIR}/${PETSC_ARCH}/include -I${PETSC_DIR}/include -I${PETSC_DIR}/include/petsc/mpiuni -DHXT_HAVE_PETSC
LDFLAGS += -L${PETSC_DIR}/${PETSC_ARCH}/lib -Wl,-rpath=${PETSC_DIR}/${PETSC_ARCH}/lib -lpetsc
else
FLUID_C+=${FLUID_C} src/linear_system.c
endif
SCONTACT_C=src/quadtree.c src/scontact.c
SCONTACT_H=src/quadtree.h src/scontact.h
libmbfluid.so : ${FLUID_C} ${FLUID_H}
${CC} ${FLUID_C} -o $@ ${CFLAGS} ${LDFLAGS} -DDIMENSION=3
libscontact2.so : ${SCONTACT_C} ${SCONTACT_H}
${CC} ${SCONTACT_C} -o $@ ${CFLAGS} ${LDFLAGS} -DDIMENSION=2
libscontact3.so : ${SCONTACT_C} ${SCONTACT_H}
${CC} ${SCONTACT_C} -o $@ ${CFLAGS} ${LDFLAGS} -DDIMENSION=3
clean :
${RM} *.so *.pyc -r __pycache__
class MshBaseElementType :
def __init__(self, name, tag, dimension) :
self.name = name
self.tag = tag
self.dimension = dimension
self.mshType = [None] * 11
self.mshTypeSerendip = [None] * 11
BaseType = [None]
Type = [None]
for t in [
('TYPE_PNT', 1, 0),
('TYPE_LIN', 2, 1),
('TYPE_TRI', 3, 2),
('TYPE_QUA', 4, 2),
('TYPE_TET', 5, 3),
('TYPE_PYR', 6, 3),
('TYPE_PRI', 7, 3),
('TYPE_HEX', 8, 3),
('TYPE_POLYG', 9, 2),
('TYPE_POLYH', 10, 3)
] :
BaseType.append(MshBaseElementType(*t))
globals()[t[0]] = BaseType[-1]
class MshElementType :
def __init__(self, name, tag, numVertices, baseType, order, serendip) :
self.name = name
self.tag = tag
self.numVertices = numVertices
self.baseType = baseType
self.order = order
self.serendip = serendip
for i, t in enumerate([
'MSH_LIN_2', 'MSH_TRI_3', 'MSH_QUA_4', 'MSH_TET_4', 'MSH_HEX_8',
'MSH_PRI_6', 'MSH_PYR_5', 'MSH_LIN_3', 'MSH_TRI_6', 'MSH_QUA_9',
'MSH_TET_10', 'MSH_HEX_27', 'MSH_PRI_18', 'MSH_PYR_14', 'MSH_PNT',
'MSH_QUA_8', 'MSH_HEX_20', 'MSH_PRI_15', 'MSH_PYR_13', 'MSH_TRI_9',
'MSH_TRI_10', 'MSH_TRI_12', 'MSH_TRI_15', 'MSH_TRI_15I', 'MSH_TRI_21',
'MSH_LIN_4', 'MSH_LIN_5', 'MSH_LIN_6', 'MSH_TET_20', 'MSH_TET_35',
'MSH_TET_56', 'MSH_TET_22', 'MSH_TET_28', 'MSH_POLYG_', 'MSH_POLYH_',
'MSH_QUA_16', 'MSH_QUA_25', 'MSH_QUA_36', 'MSH_QUA_12', 'MSH_QUA_16I',
'MSH_QUA_20', 'MSH_TRI_28', 'MSH_TRI_36', 'MSH_TRI_45', 'MSH_TRI_55',
'MSH_TRI_66', 'MSH_QUA_49', 'MSH_QUA_64', 'MSH_QUA_81', 'MSH_QUA_100',
'MSH_QUA_121', 'MSH_TRI_18', 'MSH_TRI_21I', 'MSH_TRI_24', 'MSH_TRI_27',
'MSH_TRI_30', 'MSH_QUA_24', 'MSH_QUA_28', 'MSH_QUA_32', 'MSH_QUA_36I',
'MSH_QUA_40', 'MSH_LIN_7', 'MSH_LIN_8', 'MSH_LIN_9', 'MSH_LIN_10',
'MSH_LIN_11', 'MSH_LIN_B', 'MSH_TRI_B', 'MSH_POLYG_B', 'MSH_LIN_C',
'MSH_TET_84', 'MSH_TET_120', 'MSH_TET_165', 'MSH_TET_220', 'MSH_TET_286',
None, None, None, 'MSH_TET_34', 'MSH_TET_40',
'MSH_TET_46', 'MSH_TET_52', 'MSH_TET_58', 'MSH_LIN_1', 'MSH_TRI_1',
'MSH_QUA_1', 'MSH_TET_1', 'MSH_HEX_1', 'MSH_PRI_1', 'MSH_PRI_40',
'MSH_PRI_75', 'MSH_HEX_64', 'MSH_HEX_125', 'MSH_HEX_216', 'MSH_HEX_343',
'MSH_HEX_512', 'MSH_HEX_729', 'MSH_HEX_1000', 'MSH_HEX_32', 'MSH_HEX_44',
'MSH_HEX_56', 'MSH_HEX_68', 'MSH_HEX_80', 'MSH_HEX_92', 'MSH_HEX_104',
'MSH_PRI_126', 'MSH_PRI_196', 'MSH_PRI_288', 'MSH_PRI_405', 'MSH_PRI_550',
'MSH_PRI_24', 'MSH_PRI_33', 'MSH_PRI_42', 'MSH_PRI_51', 'MSH_PRI_60',
'MSH_PRI_69', 'MSH_PRI_78', 'MSH_PYR_30', 'MSH_PYR_55', 'MSH_PYR_91',
'MSH_PYR_140', 'MSH_PYR_204', 'MSH_PYR_285', 'MSH_PYR_385', 'MSH_PYR_21',
'MSH_PYR_29', 'MSH_PYR_37', 'MSH_PYR_45', 'MSH_PYR_53', 'MSH_PYR_61',
'MSH_PYR_69', 'MSH_PYR_1', 'MSH_PNT_SUB', 'MSH_LIN_SUB', 'MSH_TRI_SUB',
'MSH_TET_SUB', 'MSH_TET_16'
]) :
if not t :
Type.append(None)
continue
w = t.split("_")
baseType = globals()["TYPE_"+w[1]]
nnodes2order = {
TYPE_LIN:(dict((o+1,o) for o in range(11)), {}),
TYPE_TRI:(dict(((o+1) * (o+2) / 2,o) for o in range(11)), dict((3 * o,o) for o in range(11))),
TYPE_QUA:(dict(((o+1) * (o+1),o) for o in range(11)), dict((4 * o,o) for o in range(11))),
TYPE_TET:(dict(((o+1) * (o+2) * (o+3) / 6,o) for o in range(11)), dict((4 + 6 * (o - 1), o) for o in range(11))),
TYPE_HEX:(dict(((o+1) * (o+1) * (o+1),o) for o in range(11)), dict((8 + 12 * (o - 1), o) for o in range(11))),
TYPE_PRI:(dict(((o+1) * (o+2) * (o+1)/2,o) for o in range(11)), dict((6 + 9 * (o - 1), o) for o in range(11))),
TYPE_PYR:(dict(((o+1) * (o+2) * (2*o +3)/6, o) for o in range(11)), dict((5 + 8 * (o - 1),o) for o in range(11)))
}
try :
numnodes = int(w[2].strip("I"))
except:
numnodes = 1 if baseType == TYPE_PNT else -1
try :
serendip = (w[2][-1] == 'I')
except :
serendip = False
if numnodes == -1 :
order = -1
elif baseType == TYPE_PNT :
order = 0
serendip = False
else:
if (not serendip) and numnodes in nnodes2order[baseType][0] :
order = nnodes2order[baseType][0][numnodes]
else :
order = nnodes2order[baseType][1][numnodes]
serendip = True
Type.append(MshElementType(t, i + 1, numnodes, baseType, order, serendip))
if serendip :
baseType.mshTypeSerendip[order] = Type[-1]
else :
baseType.mshType[order] = Type[-1]
globals()[t] = Type[-1]
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 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) :
"""
Initialize LMGC90
"""
chipy.checkDirectories()
### Set dimension in chipy for dummies ###
chipy.SetDimension(2)
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._nbDisk = chipy.DISKx_GetNbDISKx()
self._d2bMap = chipy.DISKx_GetPtrDISKx2RBDY2()
self._nbPoly = chipy.POLYG_GetNbPOLYG()
self._p2bMap = chipy.POLYG_GetPtrPOLYG2RBDY2()
self._position = np.zeros([self._nbDisk,3], 'd')
self._velocity = np.zeros([self._nbDisk,3], 'd')
self._externalF= np.zeros([self._nbDisk,3], 'd')
self._volume = np.zeros([self._nbDisk,1], 'd')
self._mass = np.zeros([self._nbDisk,1], 'd')
self._r = np.zeros([self._nbDisk,1], 'd')
for i in range(self._nbDisk):
self._r[i] = chipy.DISKx_GetContactorRadius(i+1)
self._mass[i] = chipy.RBDY2_GetBodyMass(int(self._d2bMap[i,0]))
self._volume = np.pi * self._r**2 * 2./3
self._tag2bnd = {}
self._tag2id = {}
for i in range(chipy.RBDY2_GetNbRBDY2()):
for t in range(chipy.RBDY2_GetNbContactor(i+1)):
tag = chipy.RBDY2_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-4,
'relax' : 1.,
'gsit1' : 200,
'gsit2' : 10
}
def setBoundaryVelocity(self, tag, v) :
self._tag2bnd[tag] = v
def write(self, odir, i, t) :
refRadius=np.amin(self._r)*2.
chipy.WriteOutDof(1)
chipy.WriteOutVlocRloc(1)
chipy.WritePostproFiles()
chipy.WriteDisplayFiles(1, refRadius)
#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[:,:2] = forces
chipy.IncrementStep()
for tag, v in self._tag2bnd.iteritems() :
for i in self._tag2id[tag] :
for j, iv in enumerate(v):
chipy.RBDY2_SetVlocyDrivenDof(i+1,j+1,iv)
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(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 r(self):
""" Return the radius of the particles """
return self._r
def position(self):
""" Get current position of contactors and return it """
for i in range(self._nbDisk):
self._position[i,:] = chipy.RBDY2_GetBodyVector('Coor_',int(self._d2bMap[i,0]))
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()
def scontactToLmgc90(filename, fric = 0.):
import scontact2
from pylmgc90 import pre
datbox_path = 'DATBOX'
if not os.path.isdir(datbox_path):
os.mkdir(datbox_path)
sc = scontact2.ParticleProblem()
sc.load(filename)
x = sc.position()
v = sc.velocity()
v[:] = 0.
r = sc.particles()[:, 0]
rho_particle = np.mean(sc.particles()[:, 1]/(np.pi * r**2))
avs = pre.avatars()
mat = pre.materials()
mod = pre.models()
svs = pre.see_tables()
tac = pre.tact_behavs()
mater = pre.material(name='STONE',materialType='RIGID',density=rho_particle)
model = pre.model(name='rigid', physics='MECAx', element='Rxx2D',dimension=2)
mat.addMaterial(mater)
mod.addModel(model)
for i in range(r.size) :
P = pre.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)
seg = sc.segments()
bndtags = set()
for i in range(seg.shape[0]) :
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)
exit()
tag += "_" * (5 - len(tag))
bndtags.add(tag)
t = (x0 - x1)/np.linalg.norm(x0 - x1)
n = np.array([-t[1], t[0]]) * 1e-4
vs = np.array([x0-2*n, x0, x1, x1 -2*n])
av = pre.rigidPolygon(model,mater,np.zeros([2]),color=tag,generation_type="full",vertices=vs)
av.imposeDrivenDof(component=[1,2,3],dofty='vlocy')
avs.addAvatar(av)
clb_fric = pre.tact_behav(name='iqsc0',law='IQS_CLB',fric=fric)
tac += clb_fric
svs += pre.see_table(CorpsCandidat ="RBDY2", candidat ="DISKx", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY2", antagoniste="DISKx", colorAntagoniste='INxxx',
behav=clb_fric, alert=r.min())
for tag in bndtags :
svs += pre.see_table(CorpsCandidat ="RBDY2", candidat ="DISKx", colorCandidat ='INxxx',
CorpsAntagoniste="RBDY2", antagoniste="POLYG", colorAntagoniste=tag,
behav=clb_fric, alert=r.min())
post = pre.postpro_commands()
solv = pre.postpro_command(name='SOLVER INFORMATIONS', step=1)
post.addCommand(solv)
# file writting
pre.writeBodies(avs,chemin=datbox_path)
pre.writeDrvDof(avs,chemin=datbox_path)
pre.writeDofIni(avs,chemin=datbox_path)
pre.writeModels(mod,chemin=datbox_path)
pre.writeBulkBehav(mat,chemin=datbox_path, gravy=[0.,0.,0.])
pre.writeTactBehav(tac,svs,chemin=datbox_path)
pre.writeVlocRlocIni(chemin=datbox_path)
pre.writePostpro(commands=post, parts=avs, path=datbox_path)
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()
chipy.PRPRx_UseCpCundallDetection(10)
##### 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_GetPtrSPHER2BDYTY()
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')
self._rmin=1e20
for i in range(self._nbSpher):
self._rmin = min(self._rmin,chipy.SPHER_GetContactorRadius(i+1))
self._volume[i] = np.pi * chipy.SPHER_GetContactorRadius(i+1)**3 * 4./3
self._mass[i] = chipy.RBDY3_GetMass(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, self._rmin*2.)
#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._nbSpher):
self._position[i,:] = chipy.SPHER_GetContactorCoor(i+1)[:3]
return self._position
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]))
return self._velocity[:,:3]
def externalForces(self):
""" Get an external forces array
"""
return self._externalF[:,:3]
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)