Commit 6e50a531 authored by Nathan Coppin (nathan.coppin@uclouvain.be)'s avatar Nathan Coppin (nathan.coppin@uclouvain.be)
Browse files

message

parent 7c6b89b4
Pipeline #9103 passed with stages
in 2 minutes and 33 seconds
......@@ -28,11 +28,14 @@
from migflow import fluid
from migflow import scontact
from migflow import time_integration
import numpy as np
import os
import time
import shutil
import random
import sys
numDep = int(sys.argv[1])
def genInitialPosition(filename, r, H, ly, lx, lz, rin, rhop) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure
......@@ -61,30 +64,30 @@ def genInitialPosition(filename, r, H, ly, lx, lz, rin, rhop) :
p.add_particle((x[i], y[j], z[k]), r1, (4/3)*r1**3 * np.pi * rhop,"Sand")
p.write_vtk(filename,0,0)
# Define output directory
outputdir = "output"
outputdir = str(numDep)
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
# Physical parameters
g = np.array([0,-9.81,0]) #gravity
r = 8e-3 #particles radius
r = 3e-3 #particles radius
rin = 0.01905 #central cylinder radius
rhop = 2500 #particles density
# Numerical parameters
outf = 50 #number of iterations between output files
outf = 5000 #number of iterations between output files
dt = 1e-3 #time step
tEnd = 10 #final time
tEnd = 5 + 5*dt #final time
# Geometrical parameters
ly = 0.7 #particles area height
lx = 0.15 #particles area widht
ly = 0.7 #particles area height
lx = 0.15 #particles area widht
lz = 0.05 #particles area depth
H = 1.3 #domain height
H = 1.4 #domain height
#
# PARTICLE PROBLEM
#
# Initialise particles
genInitialPosition(outputdir, r, H, ly, lx, lz, rin, rhop)
p = scontact.ParticleProblem(3,True,True)
p.set_friction_coefficient(0.3,"Sand","Sand")
p.set_friction_coefficient(0.15,"Sand","Sand")
p.set_friction_coefficient(0.3,"Sand","Steel")
p.read_vtk(outputdir,0)
# Initial time and iteration
......
import numpy as np
from migflow import scontact
from migflow import VTK
import os
directory = "depot/output3"
#Step between points where the tensor is computed
step = 0.015
#Computation radius
radius = 0.015
#Box dimensions
p = scontact.ParticleProblem(3,True,True)
p.read_vtk(directory,0)
x = p.disks()["x"]
xmin,ymin,zmin = np.min(x,axis=0)
xmax,ymax,zmax = np.max(x,axis=0)
z = np.arange(zmin,zmax,step)
x = np.arange(xmin,xmax,step)
y = np.arange(ymin,ymax,step)
x,y,z = np.meshgrid(x,y,z)
points = np.column_stack([x.flat,y.flat,z.flat])
p = scontact.ParticleProblem(3,True,True)
for filename,it,t in VTK.read_index(directory+"/particles.pvd") :
print(it)
p.read_vtk(directory,it)
stress = p.computeStressTensor(points,radius)
tensor = stress[:,[0,4,8,1,5,2]]
lmax = np.amax(abs(np.linalg.eigvalsh(stress.reshape((-1,3,3)))),axis=1)
types = np.full([points.shape[0]],1)
location = np.arange(0,points.shape[0])
cells = np.zeros((points.shape[0],1))
cells[:,0] = np.arange(0,points.shape[0],1)
VTK.write(directory+"/stress",it,t,
(types,cells.reshape([-1]),location),points,
data=(("stress",tensor.reshape([-1,6])),("lmax",lmax.reshape([-1,1]))))
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