Commit dc0c7515 authored by Hugo Ginestet's avatar Hugo Ginestet
Browse files

missing testcase

parent 015397a0
Pipeline #5452 passed with stage
in 26 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_contour_rectangle2(fic):
fic.write('Physical Surface("Domain") = {1};\n')
fic.write('Physical Line("Bottom") = {1};\n')
fic.write('Physical Line("Top") = {3};\n')
fic.write('Physical Line("Lateral") = {2,4};\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))
def creation_maillage2(name,para):
fic = open(name+'.geo',"w")
initialisation(fic,para[2],para[3])
creation_box(fic,para[0],para[1])
creation_contour_rectangle2(fic)
fic.close()
os.system("gmsh -2 {}.geo".format(name))
#os.system("rm {}.geo".format(name))
creation_maillage2('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
'''
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 and norm0>=r:
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
elif norm0 < r :
for i in range(2):
new_pos[i] += depl[i]*(norm0 - norm_max)/(r - norm_max)
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 = 10 # number of iterations between output files
dt = 100*0.5e-3 # time step
tEnd = 500
# final time
nbtot = int(tEnd/dt)
#PARTICLE PROBLEM
#
#p = scontact.ParticleProblem(2)
#p.read_vtk(outputdir,0)
def grains_v_x(x) :
return cl_vitesse[0]
def grains_v_y(x) :
return cl_vitesse[1]
#
# 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")
#p.set_boundary_velocity("Grains", cl_vitesse)
#fluid.set_open_boundary("Grains", velocity=[grains_v_x,grains_v_y])
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 =geo_ligne(1,rayon,[1,0])+geo_cercle(nb2,rayon)
#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.]
depl[0] += liste_depl[1][0]
depl[1] += liste_depl[1][1]
for k in range(2,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("***************" , np.max(np.abs(fluid.velocity_mesh())))
print(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)
t+=dt
fluid.export_vtk(outputdir,t, (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
for k in range(n):
fluid.implicit_euler(dt,newton_max_it=10)
t+=dt
fluid.export_vtk(outputdir, t, (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))
'''
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