Commit 4721733a authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

adaptation1-jon

parent 50620881
Pipeline #5438 passed with stage
in 27 seconds
from migflow import fluid
from migflow import scontact
import copy
import numpy as np
import math
import os
import time
import shutil
import random
L = 10.
lc = 0.2
h_min =0.25
h_max= 0.25
r = 1.0
para = [L,lc,h_min,h_max,[0.,0.,r]]
###################################################################################################################
"""
Création de la géométrie pour la création du maillage
"""
def creation_point(fic,Point):
fic.write('Point('+str(Point[0])+') = {'+str(Point[1])+', '+str(Point[2])+', 0.0, '+str(Point[3])+'};\n')
def creation_ligne(fic,i,Point1,Point2):
fic.write('Line('+str(i)+') = {'+str(Point1)+', '+str(Point2)+'};\n')
def creation_loop(fic,i,Ligne1,Ligne2,Ligne3,Ligne4):
fic.write('Line Loop('+str(i)+') = {'+str(Ligne1)+', '+str(Ligne2)+', '+str(Ligne3)+', ' +str(Ligne4)+'};\n\n')
def initialisation(fic,h_min,h_max):
fic.write('SetFactory("OpenCASCADE");\n')
fic.write('Mesh.CharacteristicLengthMin ='+ str(h_min) +' ;\n')
fic.write('Mesh.CharacteristicLengthMax ='+ str(h_max) +' ;\n\n')
fic.write('Mesh.MinimumCirclePoints = 1 ;\n\n')
fic.write('Mesh.MinimumCurvePoints = 2 ;\n\n')
def creation_box(fic,L,lc):
Point = [[1,-L/2.,-L/2.,lc],[2,L/2.,-L/2.,lc],[3,L/2.,L/2.,lc],[4,-L/2.,L/2.,lc]]
creation_point(fic,Point[0])
for k in range(3):
creation_point(fic,Point[k+1])
creation_ligne(fic,k+1,Point[k][0],Point[k+1][0])
creation_ligne(fic,4,Point[-1][0],Point[0][0])
creation_loop(fic,1,Point[0][0],Point[1][0],Point[2][0],Point[3][0])
fic.write('Surface(1) = {1}; \n\n')
def creation_disk(fic,i,Point):
fic.write('Disk('+str(i)+') = {'+str(Point[0])+', '+str(Point[1])+', 0., '+str(Point[2])+'};\n\n')
def creation_surface(fic,nb):
fic.write('BooleanDifference{Surface{1};Delete;}{Surface{2};Delete;}\n\n')
def creation_contour_rectangle(fic):
fic.write('Physical Surface("Domain") = {1};\n')
fic.write('Physical Line("Bottom") = {6};\n')
fic.write('Physical Line("Top") = {9};\n')
fic.write('Physical Line("Lateral") = {7,8};\n')
fic.write('Physical Line("Grains") = {5};\n\n')
fic.write('Mesh.MshFileVersion = 2;')
def creation_maillage(name,para):
fic = open(name+'.geo',"w")
initialisation(fic,para[2],para[3])
creation_box(fic,para[0],para[1])
creation_disk(fic,2,para[4])
creation_surface(fic,1)
creation_contour_rectangle(fic)
fic.close()
os.system("gmsh -2 {}.geo".format(name))
#os.system("rm {}.geo".format(name))
creation_maillage('mesh',para)
para = [L,lc,h_min,h_max]
def depl_lineaire(pos,para,depl,pos_c,dt):
l = para[0]/2
#norm0 = math.sqrt((pos[0]-pos_c[0])**2 + (pos[1]-pos_c[1])**2)
velocity = [0, 0]
'''
if abs(norm0 - pos_c[2]) < 1e-4:
for i in range(2):
pos[i] += depl[i]
velocity[i] = depl[i]/dt
'''
for i in range(2):
if pos[i] <(pos_c[i] - pos_c[2]):
pos[i] += depl[i]*(pos[i]+l)/(pos_c[i]- pos_c[2]+l)
velocity[i] = (depl[i]*(pos[i]+l)/(pos_c[i]- pos_c[2]+l))/dt
elif pos[i] >(pos_c[i] + pos_c[2]):
pos[i] += depl[i]*(pos[i]-l)/(pos_c[i]+ pos_c[2]-l)
velocity[i] = (depl[i]*(pos[i]-l)/(pos_c[i]+ pos_c[2]-l))/dt
else:
pos[i]+= depl[i]
velocity[i] = depl[i]/dt
return velocity
def depl_lineaire2(pos,para,depl,pos_c):
l = para[0]/2
norm0 = math.sqrt((pos[0])**2 + (pos[1])**2)
if abs(pos[0]-pos_c[0]) > 1e-5:
x = (pos[1]-pos_c[1])/(pos[0]-pos_c[0])
pos0 = [norm0/math.sqrt(1+x**2),norm0*x/math.sqrt(1+x**2)]
else:
pos0 = [pos[0]-pos_c[0],pos[1]-pos_c[1]]
for i in range(2):
if pos[i] <pos_c[i]:
pos[i] += depl[i]*(pos0[i]+l)/(-pos_c[2]+l)
else:
pos[i] += depl[i]*(pos0[i]-l)/(pos_c[2]-l)
return
def calcul_norm(para,depl,pos_c):
l = para[0]/2
for k in range(2):
if pos_c[0] > 0:
if pos_c[1]>0:
return min(abs(pos_c[0]-l),abs(pos_c[1]-l))
else:
return min(abs(pos_c[0]-l),abs(pos_c[1]+l))
else:
if pos_c[1]>0:
return min(abs(pos_c[0]+l),abs(pos_c[1]-l))
else:
return min(abs(pos_c[0]+l),abs(pos_c[1]+l))
def depl_radial(pos,depl,pos_c,norm_max,para,dt):
r = pos_c[2]
norm0 = math.sqrt((pos[0]-pos_c[0])**2+(pos[1]-pos_c[1])**2)
if norm0 < norm_max:
for i in range(2):
pos[i] += depl[i]*(norm0 - norm_max)/(r - norm_max)
return [depl[0]*(norm0 - norm_max)/((r - norm_max)*dt), depl[1]*(norm0 - norm_max)/((r - norm_max)*dt)]
def depl_radial0(pos0,pos,depl,pos_c,norm_max,dt):
r = pos_c[2]
norm0 = math.sqrt((pos0[0]-pos_c[0])**2+(pos0[1]-pos_c[1])**2)
new_pos = copy.deepcopy(pos0)
velocity = [0,0]
if norm0 < norm_max:
for i in range(2):
new_pos[i] += depl[i]*(norm0 - norm_max)/(r - norm_max)
if abs(r-norm0)>1e-10:
velocity[i] = (new_pos[i] - pos[i])/dt
else:
for i in range(2):
new_pos[i] = pos[i]
return [velocity,new_pos]
def geo_ligne(nb,l,vec):
L = [[0,0]]
pas = l/nb
norm =math.sqrt(vec[0]**2+vec[1]**2)
for k in range(nb):
L.append([pas*vec[0]/norm,pas*vec[1]/norm])
return L
def geo_cercle(nb,r):
theta = 2*math.pi/nb
L = [[0,0]]
for k in range(1,nb):
toto = [r*math.cos(k*theta) - r*math.cos((k-1)*theta), r*math.sin(k*theta) -r*math.sin((k-1)*theta)]
L.append(toto)
return L
def retourne_zero(pos_c,nb):
return [[(0.-pos_c[0])/nb,(0.-pos_c[1])/nb] for i in range(nb)]
#########################################################################################
outputdir = "output0"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = -9.81 # gravity
# particles radius
rhop = 1500 # particles density
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
# Numerical parameters
outf = 2 # number of iterations between output files
dt = 0.5e-2 # time step
tEnd = 5
# final time
nbtot = int(tEnd/dt)
#PARTICLE PROBLEM
#
def genInitialPosition(filename,r,L,H,para,rhop):
# Particles structure builder
p = scontact.ParticleProblem(2)
# Load mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom","Grains"])
p.add_particle((5., 5.), r, r**2 * np.pi * rhop)
p.write_vtk(filename,0,0)
#p = scontact.ParticleProblem(2)
#p.read_vtk(outputdir,0)
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
#fluid.set_wall_boundary("Grains", velocity = [10., 0.])
fluid.set_open_boundary("Grains",velocity=[10,0.])
#fluid.set_strong_boundary("Grains",0,10.)
#fluid.set_strong_boundary("Grains",1,0.)
fluid.set_wall_boundary("Top",pressure=0)
fluid.export_vtk(outputdir,0,0)
tic = time.time()
rayon = 2.
nb1 = int(nbtot/(2*(np.pi+2)))
nb2 = nbtot -2*nb1
print((nb1,nb2))
liste_depl = geo_ligne(nb1,rayon,[1,0]) + 1*geo_cercle(nb2,rayon) + geo_ligne(nb1,rayon,[-1,0])
#liste_depl = 1*(geo_ligne(nb2,rayon,[0,-5])+geo_ligne(200,rayon,[2,1]))
#liste_depl = geo_ligne(nb1,rayon,[0,-1])
n = len(liste_depl)
N = len(fluid.coordinates())
pos_centre = [0.,0.,r]
depl = copy.deepcopy(pos_centre)
#
# COMPUTATION LOOP
#
k = 0
t = 0
ii = 0
'''
pos0= copy.deepcopy(fluid.coordinates())
cl_vitesse = [0.,0.]
for k in range(n):
#norm_max = calcul_norm(para,liste_depl[k],pos_centre)
for j in range(N):
#[poubelle,fluid.coordinates()[j]] = depl_radial0(pos0[j],fluid.coordinates()[j],depl,pos_centre,L/2,dt)
[fluid.velocity_mesh()[j],fluid.coordinates()[j]] = depl_radial0(pos0[j],fluid.coordinates()[j],depl,pos_centre,L/2,dt)
cl_vitesse[0] = liste_depl[k][0]/dt
cl_vitesse[1] = liste_depl[k][1]/dt
print(cl_vitesse)
#fluid.set_open_boundary("Grains", velocity=cl_vitesse)
fluid.set_strong_boundary("Grains",0,cl_vitesse[0])
fluid.set_strong_boundary("Grains",1,cl_vitesse[1])
fluid.implicit_euler(dt,newton_max_it=10)
fluid.export_vtk(outputdir, (k+1), (k+1))
depl[0] += liste_depl[k][0]
depl[1] += liste_depl[k][1]
print('{}\n'.format(k))
cl_vitesse[0] = 0
cl_vitesse[1] = 0
fluid.set_open_boundary("Grains", velocity=cl_vitesse)
'''
for k in range(n):
fluid.implicit_euler(dt,newton_max_it=10)
fluid.export_vtk(outputdir, (n+k+1), (n+k+1))
print('{}\n'.format(k+n))
'''
liste_depl = retourne_zero(pos_centre,30)
nb = len(liste_depl)
for k in range(nb):
norm_max = calcul_norm(para,liste_depl[k],pos_centre)
for j in range(N):
#fluid.velocity_mesh()[j] = depl_radial(fluid.coordinates()[j],liste_depl[k],pos_centre,norm_max,para,dt)
fluid.velocity_mesh()[j] = depl_lineaire(fluid.coordinates()[j],para,liste_depl[k],pos_centre,dt)
#depl_lineaire2(fluid.coordinates()[j],para,liste_depl[k],pos_centre)
#fluid.implicit_euler(dt)
fluid.export_vtk(outputdir, (n+k+1), (n+k+1))
pos_centre[0] += liste_depl[k][0]
pos_centre[1] += liste_depl[k][1]
print('{}\n'.format(k))
os.system("paraview {}/fluid.pvd".format(outputdir))
'''
Supports Markdown
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