Commit 5cc7ccf3 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

sablier-2d : lmgc

parent e5c32b49
......@@ -98,20 +98,20 @@ class LmgcInterface(object):
self._freq_detect = 1
self._solver_params = { 'type' :'Stored_Delassus_Loops ',
'norm' : 'Quad ',
'conv' : 1e-5,
'conv' : 1e-4,
'relax' : 1.,
'gsit1' : 200,
'gsit2' : 200
'gsit2' : 10
}
def setBoundaryVelocity(self, tag, v) :
self._tag2bnd[tag] = v
def write(self, odir, i, t) :
def write(self, odir, i, t, refRadius=0.01) :
chipy.WriteOutDof(1)
chipy.WriteOutVlocRloc(1)
chipy.WritePostproFiles()
chipy.WriteDisplayFiles(1, 0.01)
chipy.WriteDisplayFiles(1, refRadius)
......@@ -188,7 +188,7 @@ class LmgcInterface(object):
chipy.CloseDisplayFiles()
chipy.Finalize()
def scontactToLmgc90(filename):
def scontactToLmgc90(filename, fric = 0):
import scontact2
from pylmgc90 import pre_lmgc
datbox_path = 'DATBOX'
......@@ -198,8 +198,6 @@ def scontactToLmgc90(filename):
sc = scontact2.ParticleProblem()
sc.load(filename)
fric = 0
x = sc.position()
v = sc.velocity()
v[:] = 0.
......
......@@ -18,7 +18,6 @@ if (FLTK_LIBRARIES)
endif(ENABLE_LIBPNG)
target_link_libraries(scontactplot ${EXT_LIB})
install(TARGETS scontactplot DESTINATION bin)
message("coucou ${COMP_FLAGS}")
if(COMP_FLAGS)
set_target_properties(scontactplot PROPERTIES COMPILE_FLAGS ${COMP_FLAGS})
endif(COMP_FLAGS)
......
#!/usr/bin/env python
use_lmgc=True
import pyfefp
from pyfefp.legendre import *
from pyfefp import scontact2Interface
if use_lmgc :
from pyfefp import lmgc2Interface
import scontact2
import numpy as np
from pyfefp import mesh
......@@ -16,7 +19,7 @@ def genInitialPosition(filename, r, ox, oy, lx, ly, rhop) :
y = np.arange(oy+r, oy+ly, 2*r)
for i in range(x.shape[0]) :
for j in range(y.shape[0]) :
rr=r*(0.8+random.random()*0.2)
rr=r*(0.95+random.random()*0.05)
p.addParticle((x[i], y[j]), rr, rr**2 * np.pi * rhop);
p.write(filename)
......@@ -47,11 +50,16 @@ epsilon = 1e-5 #2e-2*c*h # ?? not sure ??
print("c = %g, eps = %g" % (c, epsilon))
y0part = 0.051
y0part = 0.071
lypart = 0.025
lxpart = 0.015
genInitialPosition(outputdir + "/part-00000", r, 0., y0part, lxpart, lypart, rhop)
p = scontact2Interface.ScontactInterface(outputdir+"/part-00000")
if use_lmgc :
friction=0.1
lmgc2Interface.scontactToLmgc90(outputdir+"/part-00000", friction)
p = lmgc2Interface.LmgcInterface()
else :
p = scontact2Interface.ScontactInterface(outputdir+"/part-00000")
deltamass = p.volume().sum()*(rhop-rho)
deltapinit = deltamass*-g[1]/lxpart
......@@ -69,17 +77,17 @@ fluid.add_boundary_condition("Box", 1, 0.)
fluid.add_boundary_condition("Box", 0, 0.)
fluid.complete()
y = fluid.mesh_position()[...,1].transpose()
pinit = (y0part+lypart-y)*deltapinit/lypart
pinit[pinit<0] = 0
pinit[pinit>deltapinit] = deltapinit
fluid._dof.setField(2, pinit, fluid.solution())
#y = fluid.mesh_position()[...,1].transpose()
#pinit = (y0part+lypart-y)*deltapinit/lypart
#pinit[pinit<0] = 0
#pinit[pinit>deltapinit] = deltapinit
#fluid._dof.setField(2, pinit, fluid.solution())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
tic = time.clock()
while t < tEnd :
forces = fluid.solve() + g * p.mass()
NS = 1
NS = 5
if ii%NS == 0 :
p.iterate(dt*NS, forces)
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
......
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