yade-users team mailing list archive

[Tanvir Hossain <question556258@xxxxxxxxxxxxxxxxxxxxx>]

Question #556258 on Yade changed:
https://answers.launchpad.net/yade/+question/556258

Tanvir Hossain posted a new comment:
Hi Jan, here is the script. 
#New version, amelioration wrt 20150715 : 
#- IMPORTANT Modify the length of the vector evaluated from ndimz to ndimz. Works with Yade6 and is not adapted to the other cases
#- Add a limiter to avoid negative tau transmitted to the fluid (unphysical and source of numerical instabilities
#- Add a new configuration for the turbulentFluctuation function. Correspond to a new formulation of the C++ HydroForceEngine. Possibility to use a new turbulent fluctuation model where each particle fluctuation lifetime depends on its wall-normal position. 
#- Add an option for fluidProfile imposed case.  
#- Register the averaged velocity of the part1 and 2 in twoSize mode (Yade6)

#If activated, adapt the script to the associated computer path
Cluster = 0	
philippeComputer = 0
julienComputer = 0
lamiaComputer = 0
raphaelComputer = 0

#Import libraries
from yade import pack, plot
import math
import random as rand
import numpy as np
np.set_printoptions(threshold=numpy.nan)
import yade.timing
import pylab as pyl
from matplotlib import pyplot as pyp
import matplotlib.gridspec as gridspec
#Register the path for later use (tests + save)
scriptPath = os.path.abspath(os.path.dirname(sys.argv[-1]))


##
## Main parameters of the simulation
##
diameterPart = 6e-3	#Diameter of the particles, in meter
densPart = 2500		#density of the particles, in kg/m3
waterDepth = 7.		#Water depth in diameter
ndimz = 525  	#Number of cells in the height
from nsmp1d_yade3D525Phi061 import * 


lengthCell = 15	#Streamwise length of the periodic cell 
widthCell = 15	#Spanwise length of the periodic cell 
Nlayer = 10.		#nb of layer of particle, in diameter

fluidHeight = (Nlayer+waterDepth)*diameterPart  #Height of the flow from
the bottom of the sample


quasi2D = 0
twoSize = 0 

#If twoSize option activated, overwrite the previous input #tanvir
if twoSize:
	Nlayer1 = 10		#nb of layer of particle type 1
	Nlayer2 = 1 		#nb of layer of particle type 2
	diameterPart1 = 6e-3	#Diameter of the particles type 1, in meter
	diameterPart2 = 4e-3	#Diameter of the particles type 2, in meter
	fluidHeight = Nlayer1*diameterPart1 + Nlayer2*diameterPart2 + 7*diameterPart1	#free surface position
	diameterPart = diameterPart1 #Size 1 taken as reference

	
#Fluid library importation and precisions
phiPartMax = 0.61		#Value of the maximum volume fraction used for the thresholding


##
## Option of the simulation
##


#Imposed fluid profile? If yes, require to add a file fluid.py in the folder with inside the value of vxFluid vector + turbStress vector if you want to add turbulent fluctuations also.
imposedFluidPro = 0

#Temporary imposed profile: impose a fluid profile for 50 seconds (considered as equilibrium) and switch to the fluid resolution
initImposedFluidPro = 0


#Activation and choice of the fluid turbulent fluid velocity fluctuation. Cannot choose both. If none, will not apply fluct.
turbModel1 = 1
turbModel2 = 0

#Stress averaging: if activated, compute and save the stress tensor depth profile at each execution of measure()
stressSave = 0
pythonStressCalculation = 0	#use the python function instead of C++

#Option to put lateral walls, to break the biperiodicity and study the influence of walls. 
lateralWalls = 0

#Option for time averaging of the DEM data before fluid resolution
timeAveraging = 0 #1 activated, 0 not activated
nbAv = 5. #number of averaging between the fluid resolution

#Averaging parameter : if activated, evaluate the average using python function partConcVelocity() instead of the C++ one
pythonAverage = 0

#debug stress mode, record all the quantities necessary to study the momentum balance in detail. 
debugStress = 0


#Video
video = 1	#1, make a movie 
if video: 
	snapshotEngine = 0	#To make the movie later using mencoder (require the qt.viewer open, obtain what you see)
	paraview = 1	#To make the movie later using paraview (do not require the qt.viewer open, faster)
	fps = 1.	#frame per second #tanvir

##
## Secondary parameters of the simulation (less changed)
##

#Particles parameters
restitCoef = 0.5	#Restitution coefficient of the particles
partFrictAngle = atan(0.4)	#friction angle of the particles, in radian
slope = 0.1		#Angle of the slope in radian #tanvir

# Fluid parameters
expoRZ = 3.1		# Richardson Zaki exponent for the hindrance function of the drag
dpdx = 0.0		#pressure gradient along streamwise direction
densFluid = 1000.	#Density of the fluid
kinematicViscoFluid = 1e-6	#kinematic viscosity of the fluid

# Fluid simulation parameters (RANS), see fluidResolution function
dsig = np.zeros(ndimz)
dz =  fluidHeight/(1.0*(ndimz-1))	
sig=pyl.linspace(0e0,1e0,ndimz)
for i in range(ndimz-1):	# calcul du pas d'espace
    dsig[i] = sig[i+1] - sig[i]

#In order to fasten the script, some correspondance are supposed. Should be respected otherwise it does not work well
dtFluid = 1e-5 	#Time step for the fluid resolution
fluidResolPeriod = 1e-2	#1/fluidResolPeriod = frequency of the calculation of the fluid profile 
measurePeriod = (1e-1)*5	# We perform measurement every measurePeriod seconds #tanvir
savePeriod = 1e1 	# We save the simulation (to be able to reload) every savePeriod seconds.


#Quasi 2D option
if quasi2D ==1:
	print '\nQuasi 2D option activated\n'
	lateralWalls = 1
	widthCell = 6.5/6.
	phiPartMax = 0.51
	print 'phiPartMax = 0.51\n'


#
## Tests on the script
#
#Width check and Value of the fluid library and size of the system
if widthCell<2 and phiPartMax>0.6:
	if widthCell<=1:
		print 'Impossible to run simulation with width exactly equal to the diameter particle size !\n'
		exit()
	else:
		print '\nCareful !! Uncoherence between the 2D character and the value of phiPartMax !! \n'
#Value of ndimz and dz
if ndimz != int(round(30*fluidHeight/diameterPart)):
	print('\n!!Problem dz = d/{0} !!\nUse ndimz = {1} instead to have dz = d/30 (and change accordingly the fluid library import)\n'.format(diameterPart/fluidHeight*ndimz,int(round(30*fluidHeight/diameterPart))))
#	if quasi2D!=1:
#		exit()

#Option compatibility
if measurePeriod<fluidResolPeriod:
	print '\nNot possible to have measurePeriod<fluidResolPeriod in the present configuration, modify the script to make measurement inside measure() instead of fluidResolution() !\n'
	exit()

turbModel = 0
if turbModel1 ==1 or turbModel2==1:
	turbModel = 1
	if turbModel1 ==1 and turbModel2==1:
		print "It is not possible to have the two fluctuations model in the same time ! \n Modify turbModel1 or turbModel2"
		exit()


#Reloading simulation
reloadingSimulation = 0
if os.path.exists(scriptPath +'/paramSave.py'):
	reloadingSimulation = 0#1 Out for now, useless, it is not working

##
## Initialization of the parameters
##

# Averaging/Save
vxPart = np.zeros(ndimz) # streamwise particle velocity depth profile
vyPart = np.zeros(ndimz) # spanmwise particle velocity depth profile
vzPart = np.zeros(ndimz) # normal particle velocity depth profile
qsMean = 0		#Mean bead flux
dragTerm = np.zeros(ndimz) # drag term transmitted to the fluid vector
averageDrag = np.zeros(ndimz) # drag field vector
stressProfile = getStressProfile(1.0,ndimz,dz,0.0,np.zeros(ndimz),np.zeros(ndimz),np.zeros(ndimz))	#Matrix full of zero
positionPartX = np.zeros(len(O.bodies))	#Defined global in order to use
positionPartZ = np.zeros(len(O.bodies))	#Defined global in order to use
velocityPartX = np.zeros(len(O.bodies))	# the saving function
velocityPartZ = np.zeros(len(O.bodies))	# the saving function
if twoSize:
	phiPart1 = np.zeros(ndimz) # Concentration vector
	vxPart1 = np.zeros(ndimz) # streamwise particle velocity depth profile
	vyPart1 = np.zeros(ndimz) # spanwise particle velocity depth profile
	vzPart1 = np.zeros(ndimz) # wall-normal particle velocity depth profile
	phiPart2 = np.zeros(ndimz) # Concentration vector
	vxPart2 = np.zeros(ndimz) # streamwise particle velocity depth profile
	vyPart2 = np.zeros(ndimz) # spanwise particle velocity depth profile
	vzPart2 = np.zeros(ndimz) # wall-normal particle velocity depth profile


#Initialization of the variables necessary for timeAveraging option
if timeAveraging:
	n_dt = 1.
	vxPart_dt = np.zeros(ndimz) 
	phiPart_dt = np.zeros(ndimz)
	averageDrag_dt = np.zeros(ndimz)
	if twoSize:
		phiPart1_dt = np.zeros(ndimz)
		averageDrag1_dt = np.zeros(ndimz)
		vxPart1_dt = np.zeros(ndimz) 
		vyPart1_dt = np.zeros(ndimz)
		vzPart1_dt = np.zeros(ndimz) 
		phiPart2_dt = np.zeros(ndimz)
		averageDrag2_dt = np.zeros(ndimz)
		vxPart2_dt = np.zeros(ndimz) 
		vyPart2_dt = np.zeros(ndimz)
		vzPart2_dt = np.zeros(ndimz)


#Parameter initialized only if it is the initial simulation, otherwise they are charged
if reloadingSimulation !=1:
	#Parameters which should be loaded and not initiallized if we recover the simulation
	phiPart = np.zeros(ndimz) # Concentration vector
	turbulentViscosity = np.zeros(ndimz)	# Turbulent viscocity
	waterDepth = 0.0	#Water depth : height of the "clear" water
	vxFluid = np.zeros(ndimz)	# Fluid velocity at time step n
	turbStress = np.zeros(ndimz)	#Turbulent stress vector
	fileNumber = 0	#counter to save the files


###########################
## PARAMETERS DECLARATION##
###########################

#Geometrical configuration, define useful quantities
height = 5*fluidHeight	#heigth of the periodic cell, bigger than the real height to take into account jumps
length = lengthCell*diameterPart 	#length of the stream
width = widthCell*diameterPart  #width of the stream 
gravityVector = Vector3(9.81*sin(slope),0.0,-9.81*cos(slope))
gravityVectorApplied = Vector3(0.0,0.0,-9.81*cos(slope))
groundPosition = height/4.0 
VLayer = dz*length*width	#Volume of a layer
partVolume = pi/6.*pow(diameterPart,3)
if twoSize:
	partVolume1 = partVolume 
	partVolume2 = pi/6.*pow(diameterPart2,3)
	
#Particles contact law parameters
maxPressure = (densPart-densFluid)*phiPartMax*Nlayer*diameterPart*abs(gravityVector[2]) #Estimated max particle pressure ("hydrostatic")
normalStiffness = maxPressure*diameterPart*1e4 #Evaluate the minimal normal stiffness to be in the rigid particle limit (cf Roux and Combe 2002)
youngMod = normalStiffness/diameterPart	#Young modulus of the particles from the stiffness wanted. Need that in order to use the restitution coefficient definition +  more practical when polydisperse sample. 
poissonRatio = 0.5	#poisson's ratio of the particles. Classical values, not much influence (see ref in Da Cruz et al 2005)
if twoSize:
	youngMod1 = normalStiffness/diameterPart1
	youngMod2 = normalStiffness/diameterPart2
	O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod1, poisson=poissonRatio, density=densPart, frictionAngle=partFrictAngle, label='Mat1'))  
	O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod2, poisson=poissonRatio, density=densPart, frictionAngle=partFrictAngle, label='Mat2'))  
#Material
O.materials.append(ViscElMat(en=restitCoef, et=0., young=youngMod, poisson=poissonRatio, density=densPart, frictionAngle=partFrictAngle, label='Mat'))  


# if it is the initial simulation create the framework and the simulation loop, otherwise load it
if reloadingSimulation !=1:
	##########################################
	## FRAMEWORK CREATION (sphere, wall,...)##
	##########################################

	#Definition of the semi-periodic cell
	O.periodic = True 
	O.cell.setBox(length,width,height)

	#To be compatible with lateral walls
	leftLimitY = 0.0
	rightLimitY = width
	centerLimitY = width/2.0
	#If lateralWalls, add the walls and increase the cell size to avoid particle touching the walls to appear on the other side
	if lateralWalls:
		#Warn the user
		print '\nlateralWalls option activated: mono-periodic boundary condition only !\n'
		#Increase the cell size
		O.cell.setBox(length,width+2*diameterPart,height)
		#Modify accordingly the position of the center of the cell and the wall right and left position
		leftLimitY = diameterPart
		rightLimitY = width+diameterPart
 		centerLimitY = diameterPart + width/2.0
		#Define the wall and add to the simulation
		sidePlaneL = box(center= (length/2.0,leftLimitY,height/2.0),extents=(2000,0,height*10),fixed=True,wire = True,color = (1,0,0), material = 'Mat')
		sidePlaneR = box(center= (length/2.0,rightLimitY,height/2.0),extents=(2000,0,height*10.0),fixed=True,wire=True, material = 'Mat',color = (0,0,1))
		O.bodies.append([sidePlaneR,sidePlaneL])


	if quasi2D == 1:
		# Regular arrangement of spheres sticked at the bottom with random height
		a = range(0,int(length/(diameterPart))) #The length is divided in diameter of particle
		for b in a:      #loop creating a set of sphere sticked at the bottom with a (uniform) random altitude comprised between 0.5 (diameter/12) and 5.5mm (11diameter/12) with steps of 0.5mm. The repartition is made around groundPosition. Exactly the same set up as the one from the experiment. 
			n =  rand.randrange(0,12,1)/12.0*diameterPart     #Define a number between 0 and 11/12 diameter with steps of 1/12 diameter (0.5mm in the experiment)      
			O.bodies.append(sphere((b*diameterPart,centerLimitY,groundPosition - 11*diameterPart/12.0/2.0 + n),diameterPart/2.,color=(0,0,0),fixed = True,material = 'Mat'))

	else:
		#Execute the file containing the database which is a 2D list containing the positions and radius of the particles composing the random fixed bed. Each position x,y (integer) contain a list of particle with position and radius at this place
		if Cluster==1:
			pathFixedBot = '/home/maurin/'
		elif philippeComputer==1:
			pathFixedBot = '/home/philippe/fixedBot/'
		elif tanvir==1:
			pathFixedBot = '/home/tanvir/output/'
		elif julienComputer==1:
			pathFixedBot = '/home/julien/Yade/scriptsRaphael/'
		elif lamiaComputer==1:
			pathFixedBot = '/home/lamia/periodicStream/'
		elif raphaelComputer:
			pathFixedBot = '/home/raphael/these/scripts/twoWayCoupling/'
		else:
			print '\nNo precised computer used ! Cannot resolve the fixed bottom file path accordingly ! Precise the computer used at line 12\n'
			exit()

		if diameterPart==6e-3:
			if lengthCell>= 10:
					execfile(pathFixedBot+'fixedBotDataLongMONO.py')
			else:
					execfile(pathFixedBot+'fixedBotDataShort.py')    
					print '\n Use the short database for random fixed bottom generation ! \n Be careful, polydisperse bottom !\n'
		elif diameterPart==3e-3:
			execfile(pathFixedBot+'fixedBotDataMONOd3.py')
		elif diameterPart==12e-3:
			execfile(pathFixedBot+'fixedBotDataMONOd12.py')
		else:
			print '\nno fixed bottom of the right size !! quit the simulation ! \n'
			exit()

		xMax = len(fixedBottomData)	# Maximum size of the data base along x axis
		yMax = len(fixedBottomData[0])	# Maximum size of the data base along y axis
		x = rand.randrange(0,xMax,1)	# Generate a uniform random number between 0 and the size of the data base along x (in diameterPart)
		y = rand.randrange(0,yMax,1)	# Generate a uniform random number between 0 and the size of the data base along y (in diameterPart)
		#Starting from the random x and y position, loop over the length and width of the require fixed bed increasing x and y position (in diameterPart). Create then a sphere starting from 0 for x and y. 
		for i in range(0,lengthCell):
			for j in range(0,int(widthCell)):
				for k in fixedBottomData[(i+x)%xMax][(j+y)%yMax]:	
					O.bodies.append(sphere(center = ((i+k[0]/diameterPart-int(k[0]/diameterPart))*diameterPart,(j+k[1]/diameterPart-int(k[1]/diameterPart))*diameterPart+leftLimitY,groundPosition+k[2]),radius = k[3],material = 'Mat',fixed = True,color = (0,0,0)))


	# Ground reference and Wall on the side
	lowPlane = box(center= (length/2.0,centerLimitY,groundPosition),extents=(200,200,0),fixed=True,wire=False,color = (0.,1.,0.),material = 'Mat')  #Build a plane to have a reference for the eyes
	WaterSurface = box(center= (length/2.0,centerLimitY,groundPosition+fluidHeight),extents=(2000,width/2.0,0),fixed=True,wire=False,color = (0,0,1),material = 'Mat',mask = 0)  #Build a plane to have a reference for the eyes of the position of the water surface. No interaction with the sphere
	O.bodies.append([lowPlane,WaterSurface]) #add to simulation

	#Count the initial number of particles (for later check)
	initNumber = len(O.bodies)

	#Create a loose cloud of particle inside the cell
	partCloud = pack.SpherePack()

	if twoSize:
		diameterPart = diameterPart1
		partVolume = partVolume1
		Nlayer = Nlayer1

	#Define the number of particle from the density of particle per unit length
	if quasi2D==1:#to keep the correspondance with previous work where I putter NbPart = Nl*L/d
		partNumber =  int(Nlayer*lengthCell)
	else:
		partNumber = int(Nlayer*phiPartMax*diameterPart*length*width/partVolume) #Volume of beads
	#Define the deposition height considering that the packing realised by make cloud is 0.2
	depositionHeight = Nlayer*phiPartMax/0.2*diameterPart #Consider that the packing realised by make cloud is 0.2
	#Create a cloud of partNumber particles
	partCloud.makeCloud(minCorner=(0,leftLimitY,groundPosition+diameterPart+1e-4),maxCorner=(length,rightLimitY,groundPosition+depositionHeight),rRelFuzz=0., rMean=diameterPart/2.0, num = partNumber)
	#Send this packing to simulation with material Mat
	partCloud.toSimulation(material='Mat')
	partNumber += initNumber 

	#Same for a second size of particles
	if twoSize:
		partCloud2 = pack.SpherePack()
		if quasi2D==1:#to keep the correspondance with previous work where I putter NbPart = Nl*L/d
			partNumber2 =  int(Nlayer2*lengthCell*diameterPart1/diameterPart2)
		else:
			partNumber2 = int(Nlayer2*phiPartMax*diameterPart2*length*width/partVolume2) #Volume of beads
		#Define the deposition height considering that the packing realised by make cloud is 0.2
		depositionHeight2 =Nlayer2*diameterPart2*phiPartMax/0.2
		partCloud2.makeCloud(minCorner=(0,leftLimitY,groundPosition+depositionHeight),maxCorner=(length,rightLimitY,groundPosition+depositionHeight+depositionHeight2),rRelFuzz=0., rMean=diameterPart2/2.0, num = partNumber2)
		partCloud2.toSimulation(material='Mat2',color = (0,0,1)) #send this packing to simulation with material Mat2
		#Increment the counters for check and deposition time input
		depositionHeight+=depositionHeight2
		partNumber+=partNumber2
	#Evaluate the deposition time considering the free-fall time of the highest particle to the ground (overestimation)
	depoTime = sqrt(depositionHeight*2/abs(gravityVector[2]))

	#Check
	if len(O.bodies)<partNumber:
		print '\nProblem : too low number of beads in the sample compared with what was asked !\n'
		exit()


	# Collect the ids of the spheres which are dynamic and to which we should add the hydroforce
	dv =1e-2	#Scale of the velocity imposed to the particle, no influence if not extreme case en = 1
	idApplyForce = []
	for b in O.bodies: 
		if isinstance(b.shape,Sphere) and b.dynamic:
			if widthCell<2:	#For the quasi 2D case, need to push slightly the spheres randomly in order to break the 2D configuration
				deltaV = - dv + rand.uniform(0,2*dv)
				b.state.vel[1]+=deltaV
			idApplyForce+=[b.id]

	#########################
	#### SIMULATION LOOP#####
	#########################

	O.engines = [
		# Reset the forces
		ForceResetter(),
		# Detect the potential contacts
		InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection',allowBiggerThanPeriod = True),
		# Calculate the different interactions
		InteractionLoop(
	   	[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
	   	[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
	   	[Law2_ScGeom_ViscElPhys_Basic()]
		,label = 'interactionLoop'),				
		#Apply an hydrodynamic force to the particles
		HydroForceEngine(densFluid = densFluid,viscoDyn = kinematicViscoFluid*densFluid,zRef = groundPosition,gravity = gravityVectorApplied,deltaZ = dz,expoRZ = expoRZ,lift = False,nCell = ndimz,vCell = length*width*dz ,vxFluid = vxFluid,phiPart = phiPart,vxPart = vxPart,ids = idApplyForce, label = 'addFluidForce', dead = True),
		#Time averaging of the DEM data for fluid resolution
		PyRunner(command = 'timeAveragingDEMdata()', virtPeriod = fluidResolPeriod/nbAv, label = 'timeAve', dead = True),
		#Calculate the fluid profile using the routine
		PyRunner(command = 'fluidResolution(fluidResolPeriod)', virtPeriod = fluidResolPeriod, label = 'fluidResol', dead = True),
		#Calculate the turbulent fluctuations from a DRW random walk model
		PyRunner(command = 'turbulentFluctuation()', virtPeriod = 0.1, label = 'turbFluct', dead = True),
		#Measurement, output files
		PyRunner(command = 'measure()', virtPeriod = measurePeriod, label = 'measurement', dead = True),
		# Check if the packing is stabilized, if yes activate the hydro force on the grains and the slope.
		PyRunner(command='gravityDeposition(depoTime)',virtPeriod = 0.01,label = 'packing'),
		#GlobalStiffnessTimeStepper, determine the time step
		GlobalStiffnessTimeStepper(defaultDt = 1e-4, viscEl = True,timestepSafetyCoefficient = 0.7,  label = 'GSTS'),
		# Integrate the equation and calculate the new position/velocities...
		NewtonIntegrator(damping=0.2, gravity=gravityVector, label='newton')
		]
	#save the initial configuration to be able to recharge the simulation starting configuration easily
	O.saveTmp()
	#run
	O.run()

	#Add engines to make a video
	if video:
		if snapshotEngine:
			print '\n video snapshot activated : takes snapshot of the qt.view as you see it\n'
			O.pause()
			O.engines += [SnapshotEngine(virtPeriod=1.0/fps,fileBase=scriptPath+'/figure/',label='snapshooter')]
		elif paraview:
			print '\n video paraview activated : register vtk file for later use of paraview (video...)\n'
			O.engines += [VTKRecorder(virtPeriod=1/fps,recorders=['spheres','colors'],fileName=scriptPath+'/Paraview/p1-')]
		else: 
			print '\nBe careful the video is not activated : need to specify snapshotEngine or paraview value.\n'
			O.pause()

	#Add a pyRunner to re-activate the fluctuations at 50s
	if initImposedFluidPro:
		O.engines += [PyRunner(command='activateFluidResol()',virtPeriod = 10,label = 'actFlResol')]

# else, it is a reloading of simulation, we load the necessary parameters and the simulation
else :
	execfile(scriptPath +'/paramSave.py')	#Execute the file containing the parameter registered previously
	O.load(scriptPath + '/save.xml')	#load the simulation where it stopped
	print '\n RELOAD THE PREVIOUS SIMULATION\n'

if imposedFluidPro==1 or initImposedFluidPro==1:
	execfile('fluid.py')	#Import the fluid profile (+turbStress if necessary)


###################################################################################################################################
###################################################  FUNCTION DEFINITION  #########################################################
###################################################################################################################################


######								           ######
### LET THE TIME FOR THE GRAVITY DEPOSITION AND ACTIVATE THE FLUID AT THE END ###
######								           ######
def gravityDeposition(lim):
	#At the very start, unbalanced force can be low as there is only few contacts, but it does not mean the packing is stable	
	#The rest will be run only if O.time>lim = estimation of the deposition time
	if O.time<lim : return
	else :		
		newton.damping = 0.0
	   	packing.dead = True	# Remove the present engine for the following

		addFluidForce.dead = False	# Activate the fluid force
		fluidResol.dead = False	# Activate the fluid profile resolution
#		O.timingEnabled = True	# Activate a function to check the time taken by the different function

		# If option activated,  Activate the turbulent fluctuation calculation
		if turbModel==1:
			turbFluct.dead = False
		# If option activated, activate the function to perform cumulated time averaging on DEM for fluidResol
		if timeAveraging:
			timeAve.dead = False
		#If twoSize C++ option activated, express the size of the two particles in the engine for later averaging
		if twoSize and pythonAverage!=1:  
			addFluidForce.radiusPart1 = diameterPart1/2.
			addFluidForce.radiusPart2 = diameterPart2/2.
			addFluidForce.twoSize = True

		if imposedFluidPro==1 or initImposedFluidPro==1:
			fluidResol.dead = True	#Force no fluid resolution
			addFluidForce.vxFluid = vxFluid #pass the fluid velocity vector to add drag force value
			addFluidForce.simplifiedReynoldStresses = turbStress	#and reynolds stress tensor for fluctuations
		else:
			#Initialization of the fluid resolution
			fluidResolution(fluidResolPeriod)

		measurement.dead = False	# Activate the Measurement
	return
###############
#########################################


	
#############								###################
#  Activate the fluid resolution. Used only when initImposedFluidPro option activated	###
#############								###################
def activateFluidResol():
	global initImposedFluidPro
	if O.time>50:
		fluidResol.dead = False	#Activate the fluid resolution
		initImposedFluidPro = 0	#Put the option to zero, to avoid confusion in the script
		actFlResol.dead = True	#Kill the engine

#######################################################



###################################################################################################
##############	Function to time-average the DEM data transmitted to fluid resolution	###########
###################################################################################################
def timeAveragingDEMdata():
	global phiPart_dt
	global vxPart_dt
	global averageDrag_dt
	global phiPart1_dt
	global averageDrag1_dt
	global vxPart1_dt
	global vyPart1_dt
	global vzPart1_dt
	global phiPart2_dt
	global averageDrag2_dt
	global vxPart2_dt
	global vyPart2_dt
	global vzPart2_dt
	global n_dt

	#Re-initialized at each new averaging step
	if n_dt==0:
		vxPart_dt = np.zeros(ndimz)
		phiPart_dt = np.zeros(ndimz)
		averageDrag_dt = np.zeros(ndimz)
		if twoSize:
			phiPart1_dt = np.zeros(ndimz)
			averageDrag1_dt = np.zeros(ndimz)
			vxPart1_dt = np.zeros(ndimz)
			vyPart1_dt = np.zeros(ndimz)
			vzPart1_dt = np.zeros(ndimz)
			phiPart2_dt = np.zeros(ndimz)
			averageDrag2_dt = np.zeros(ndimz)
			vxPart2_dt = np.zeros(ndimz)
			vyPart2_dt = np.zeros(ndimz)
			vzPart2_dt = np.zeros(ndimz)
	
	#Evaluate the average vx,vy,vz,phi,drag profiles and store it in addFluidForce
	addFluidForce.averageProfile()
	
	#Store the values in 
	phiPart_dt += np.array(addFluidForce.phiPart)
	vxPart_dt += np.array(addFluidForce.vxPart)
	averageDrag_dt += np.array(addFluidForce.averageDrag)
	if twoSize:
		phiPart1_dt += np.array(addFluidForce.phiPart1)
		averageDrag1_dt += np.array(addFluidForce.averageDrag1)
		vxPart1_dt += np.array(addFluidForce.vxPart1)
		vyPart1_dt += np.array(addFluidForce.vyPart1)
		vzPart1_dt += np.array(addFluidForce.vzPart1)
		phiPart2_dt += np.array(addFluidForce.phiPart2)
		averageDrag2_dt += np.array(addFluidForce.averageDrag2)
		vxPart2_dt += np.array(addFluidForce.vxPart2)
		vyPart2_dt += np.array(addFluidForce.vyPart2)
		vzPart2_dt += np.array(addFluidForce.vzPart2)

        n_dt+=1.

###############
#########################################


################					  #########################
########## EVALUATE THE TURBULENT FLUCTUATION ASSOCIATED TO EACH PARTICLE##########
#################					  #########################
# Argument: ndimz, turbulentViscosity, vxFluid, phiPart,
# Modification : turbFluct.virtPeriod, addFluidForce.bedElevation, addFluidForce.simplifiedReynoldStresses, addFluidForce.velFluct
addFluidForce.simplifiedReynoldStresses = np.zeros(ndimz)
addFluidForce.bedElevation = 0. 
addFluidForce.turbulentFluctuation() 
addFluidForce.fluctTime = np.zeros(len(O.bodies))
def turbulentFluctuation():
	global waterDepth
	global turbStress

	#For stability requirement at the initialization stage
	if O.time<depoTime+0.5:
		print 'No turbulent fluctuation in the initialization process for stability reasons!'
		turbFluct.virtPeriod = 0.5
	else:
		# Evaluate nBed, the position of the bed which is assumed to be located around the first maximum of concentration when considering decreasing z.
		nBed = 0.
		bedElevation = 0.
		for n in range(1,ndimz):
			# if there is a peak and its value is superior to 0.5, we consider it to be the position of the bed
			if phiPart[ndimz - n] < phiPart[ndimz - n-1] and phiPart[ndimz - n] > 0.5 :
				nBed = ndimz - n
				waterDepth = (ndimz-1 - nBed)*dz
				bedElevation = fluidHeight - waterDepth	 #Evaluate the bed elevation for the following
				break
		#(Re)Define the bed elevation over which fluid turbulent fluctuations will be applied.
		addFluidForce.bedElevation = bedElevation	

		#First turbulent fluctuation model: a unique constant lifetime for the turbulent fluctuation, flucTimeScale
		if turbModel1:
			vMeanAboveBed = sum(vxFluid[nBed:])/(ndimz-nBed)	# fluid elocity scale in the water depth
			flucTimeScale = waterDepth/vMeanAboveBed	# time scale of the fluctuation w_d/v, eddy turn over time
			# New evaluation of the random fluid velocity fluctuation for each particle. 
			addFluidForce.turbulentFluctuation() 
			turbFluct.virtPeriod = flucTimeScale	#Actualize when will be calculated the next fluctuations. 

		elif turbModel2:
			# New evaluation of the random fluid velocity fluctuation for each particle. 
			addFluidForce.turbulentFluctuationBIS() 
			turbFluct.virtPeriod = min(addFluidForce.fluctTime)#Actualize when will be calculated the next fluct. 
			addFluidForce.dtFluct = min(addFluidForce.fluctTime)#Actualize the time step of the function 

		if O.time<depoTime+1: #To avoid important virtPeriod at the start, impose it for 1s. Fasten the transient.
			turbFluct.virtPeriod = 0.1
			if turbModel2: addFluidForce.dtFluct=0.1


###############
#########################################


######					  ######
### CALL THE FLUID PROFILE RESOLUTION ROUTINE###
######					  ######
#Arguments : ndimz, vxFluid, vxPart, dragTerm, densFluid, fluidHeight, sig, dsig, diameterPart, phiPartFluid, kinematicViscoFluid, densPart, dpdx, slope, tfin,d tFluid
#Modify : vxFluid, addFluidForce.vxFluid
testFluidResol = 0	# Test variable for the fluid resolution
def fluidResolution(tfin):
	global phiPart
	global vxPart
	global vyPart
	global vzPart
	global vxFluid
	global averageDrag
	global turbulentViscosity
	global testFluidResol
	global dragTerm
	global taufsi
	global n_dt
	global turbStress
	if twoSize:
		global averageDrag1
		global phiPart1
		global vxPart1
		global vyPart1
		global vzPart1
		global averageDrag2
		global phiPart2
		global vxPart2
		global vyPart2
		global vzPart2
	
	#Two different way of evaluating the average profiles: phiPart, vxPart, averageDrag, (+ same for 1 and 2 with twoSize)
	#From python
	if pythonAverage==1 or debugStress==1:
		averageProfile()
	#Or from C++
	elif timeAveraging:#with multiple time averaging made previously in timeAveragingDEMdata
		#Normalize the time averaged data
		vxPart = vxPart_dt/n_dt
		phiPart = phiPart_dt/n_dt
		averageDrag = averageDrag_dt/n_dt
		if twoSize:
			phiPart1 = phiPart1_dt/n_dt
			averageDrag1 = averageDrag1_dt/n_dt
			vxPart1 = vxPart1_dt/n_dt
			vyPart1 = vyPart1_dt/n_dt
			vzPart1 = vzPart1_dt/n_dt
			phiPart2 = phiPart2_dt/n_dt
			averageDrag2 = averageDrag2_dt/n_dt
			vxPart2 = vxPart2_dt/n_dt
			vyPart2 = vyPart2_dt/n_dt
			vzPart2 = vzPart2_dt/n_dt

		#And report the other one (for usual measurement)
		vyPart = np.array(addFluidForce.vyPart)
		vzPart = np.array(addFluidForce.vzPart)
		n_dt = 0

	else:#with simple time averaging
		if raphaelComputer!=1:
			#Evaluate the average vx,vy,vz,phi,drag profiles and store it in addFluidForce
			addFluidForce.averageProfile()
		else:#On old ubuntu/Yade version, needs to do it differently:
			#At next time step, Evaluate the average vx,vy,vz,phi,drag profiles and store it in addFluidForce
			if testFluidResol == 0:
				addFluidForce.activateAverage = True	#Activate the average for next time step
				fluidResol.iterPeriod = 1	#For it to be re-executed at next time step
				testFluidResol = 1
				return


		## Import the necessary average variables for the fluid resolution from addFluidForce : phiPart,vxPart, averageDrag
		vxPart = np.array(addFluidForce.vxPart)
		vyPart = np.array(addFluidForce.vyPart)
		vzPart = np.array(addFluidForce.vzPart)
		phiPart = np.array(addFluidForce.phiPart)
		averageDrag = np.array(addFluidForce.averageDrag)

		#Back to the usual resolution period for next step
		fluidResol.iterPeriod = -1	#For it not to be re-executed next time
		fluidResol.virtPeriod = fluidResolPeriod	#Back to the normal period 
		testFluidResol = 0	#Re-initialize the test

		#If two particle size in the simulation, store complementary necessary average
		if twoSize:
			vCell = dz*width*length
			#C++ function evaluating the depth profiles of solid volume fraction and solid velocities
			phiPart1=np.array(addFluidForce.phiPart1)
			phiPart2=np.array(addFluidForce.phiPart2)
			averageDrag1=np.array(addFluidForce.averageDrag1)
			averageDrag2=np.array(addFluidForce.averageDrag2)
			vxPart1=np.array(addFluidForce.vxPart1)
			vxPart2=np.array(addFluidForce.vxPart2)
			vyPart1=np.array(addFluidForce.vyPart1)
			vyPart2=np.array(addFluidForce.vyPart2)
			vzPart1=np.array(addFluidForce.vzPart1)
			vzPart2=np.array(addFluidForce.vzPart2)

	#
	## Prepare the variables for the fluid resolution
	#

	# Threshold the solid vol. fraction in the bed to remove the oscillations due to layering and avoid turbulence in the bed
	phiPartFluid =  np.zeros(ndimz)
	for n in range(1,ndimz+1):
		# If the concentration inside the bed reach the maximum concentration (imposed in the fluid solver). Then all the values of concentration under this one are going to be thresholded at phiPartMax, in order to damp the turbulence inside the bed. 
		if phiPart[ndimz - n] > phiPartMax:
			for i in range(n,ndimz+1):
					phiPartFluid[ndimz - i] = phiPartMax
			break
		#Otherwise we report the normal value. 
		phiPartFluid[ndimz - n] = phiPart[ndimz - n]

	#Density of particle per unit volume N/Vcell:
	nPart = phiPart[:ndimz]/partVolume

	#Calculate the momentum transfer from the drag force interaction : N<fd>/Vcell
	dragTerm = nPart*averageDrag[:ndimz]
	if twoSize: #if two size, sum the two contributions n1<f1> + n2 <f2>, formulation necessary to keep energy conservation.
		dragTerm = phiPart1[:ndimz]/partVolume1*averageDrag1[:ndimz] + phiPart2[:ndimz]/partVolume2*averageDrag2[:ndimz]

	#Create Taufsi/rhof = dragTerm/(rhof(vf-vxp)) to transmit to the fluid code
	taufsi = np.zeros(ndimz)
	for i in range(1,ndimz):
		urel = max(abs(vxFluid[i] - vxPart[i]),1e-5)	#Limit the value to avoid division by 0
		#Limit the max value of taufsi to the fluid resolution limit, i.e. 1/(10dt) and required positive (charact. time)
		taufsi[i] = max(0,min(dragTerm[i]/urel/densFluid, pow(10*dtFluid,-1)))	
	
	#
	## Call the routine for fluid resolution
	#
	[vxFluid_np,turbulentViscosity] = nsmp1d_yade(fluidHeight,sig,dsig,diameterPart,vxFluid,1e0-phiPartFluid[0:ndimz],densFluid,kinematicViscoFluid,vxPart[0:ndimz],phiPartFluid[0:ndimz],densPart,dpdx,slope,tfin,dtFluid,taufsi)


	#Debugstress option
	if debugStress==1:
		global taufsiBIS
		taufsiBIS = np.zeros(ndimz)
		for i in range(1,ndimz):
			urel = max(abs(vxFluid[i] - vxPart[i]),1e-5)
			taufsiBIS[i] = min(dragTermBIS[i]/urel/densFluid, pow(10*dtFluid,-1))	
			global vxFluid_before
			vxFluid_before = vxFluid
			vxFluid = vxFluid_np
			path = scriptPath + '/'+ str(fileNumber)+'_debug.py'
			globalParam =  ['vxPart','vxPartBIS','phiPart','phiPartBIS','vxFluid','vxFluid_before', 'averageDrag','averageDragBIS','dragTerm','dragTermBIS','taufsi','taufsiBIS']
			Save(path, globalParam)
			

	#update
	vxFluid = vxFluid_np
	
	#Modify the velocity in the HydroForceEngine which apply the force
	addFluidForce.vxFluid = vxFluid


	#Evaluate the turbulent shear stress for the turbulent fluctuations model and Shields number evaluation
	turbStress = np.zeros(ndimz)
	n = 0
	for nu_t in turbulentViscosity[:-2]:	# Turbulent viscosity is global, calculated in fluidResolution
		turbStress[n] = nu_t*(vxFluid[n+1]-vxFluid[n])/dz #tau/rho
		n+=1
	addFluidForce.simplifiedReynoldStresses = turbStress	#Actualize the Reynolds stress tensor (/rho)value


###############
#########################################


#######		      ########
###	    OUTPUT	   ###
#######		      ########
#Function which register the different quantities asked for post processing
#Arguments : len(O.bodies), O.bodies, scriptPath, fileNumber, measurePeriod, savePeriod, qsMean, phiPart, vxPart, vxFluid, turbulentViscosity, dragTerm, waterDepth
#Modify : fileNumber, 'paramSave.py', 'save.xml'
def measure():
	global fileNumber
	global qsMean
	global stressProfile
	global vxPart
	global vyPart
	global vzPart
	global phiPart
	global averageDrag
	global phiPart1
	global phiPart2
	global averageDrag1
	global averageDrag2
	global vxPart1
	global vxPart2
	global vyPart1
	global vyPart2
	global vzPart1
	global vzPart2

	if imposedFluidPro==1 or initImposedFluidPro==1:#Then you should evaluate the average (as not done in fluidResolution)
		if pythonAverage==1 or debugStress==1:
			averageProfile()
		#Or from C++
		else:
			#Evaluate the average vx,vy,vz,phi,drag profiles and store it in addFluidForce
			addFluidForce.averageProfile()

			#Replace it in the global python variables
			vxPart = np.array(addFluidForce.vxPart)
			vyPart = np.array(addFluidForce.vyPart)
			vzPart = np.array(addFluidForce.vzPart)
			phiPart = np.array(addFluidForce.phiPart)
			averageDrag = np.array(addFluidForce.averageDrag)

			#If two particle size in the simulation, store complementary necessary average
			if twoSize:
				vCell = dz*width*length
				#C++ function evaluating the depth profiles of solid volume fraction and solid velocities
				phiPart1=np.array(addFluidForce.phiPart1)
				phiPart2=np.array(addFluidForce.phiPart2)
				averageDrag1=np.array(addFluidForce.averageDrag1)
				averageDrag2=np.array(addFluidForce.averageDrag2)
				vxPart1=np.array(addFluidForce.vxPart1)
				vxPart2=np.array(addFluidForce.vxPart2)
				vyPart1=np.array(addFluidForce.vyPart1)
				vyPart2=np.array(addFluidForce.vyPart2)
				vzPart1=np.array(addFluidForce.vzPart1)
				vzPart2=np.array(addFluidForce.vzPart2)


	# Save all the following variables/split in two file for faster post processing

	if stressSave == 1:
		if pythonStressCalculation:
			stressProfile = stressFieldCalculation(ndimz,dz*width*length,vxPart,dz)
		else:
			stressProfile = getStressProfile(dz*width*length,ndimz,dz,groundPosition,vxPart,vyPart,vzPart)
		for n in range(0,ndimz):
			stressProfile[0][n] = np.array(stressProfile[0][n]).reshape(1,9)[0]
			stressProfile[1][n] = np.array(stressProfile[1][n]).reshape(1,9)[0]
		globalParam =  ['stressProfile']
		path = scriptPath + '/'+ str(fileNumber)+'_stress.py'
		Save(path, globalParam)

	#if necessary output file with the positions and velocities of the particles
	if 0:
		global positionPartX
		global positionPartZ
		global velocityPartX
		global velocityPartZ
		positionPartX = np.zeros(len(O.bodies))	#Defined global in order to use
		positionPartZ = np.zeros(len(O.bodies))	#Defined global in order to use
		velocityPartX = np.zeros(len(O.bodies))	# the saving function
		velocityPartZ = np.zeros(len(O.bodies))	# the saving function
		for b in O.bodies:
			if isinstance(b.shape,Sphere):   # If it is a sphere
				positionPartX[b.id] =  b.state.pos[0]
				positionPartZ[b.id] =  b.state.pos[2]
				velocityPartX[b.id] =  b.state.vel[0]
				velocityPartZ[b.id] =  b.state.vel[2]
		path = scriptPath + '/'+ str(fileNumber)+'_2.py'
		globalParam =  ['positionPartX', 'positionPartZ','velocityPartX', 'velocityPartZ']
		Save(path, globalParam)

	#Evaluate the dimensionless sediment rate to save
	qsMean = sum(phiPart*vxPart)*dz/sqrt((densPart/densFluid - 1)*abs(gravityVector[2])*pow(diameterPart,3))

	path = scriptPath + '/'+ str(fileNumber)+'.py'
	globalParam =  ['qsMean','phiPart','vxPart','vyPart','vzPart','vxFluid','turbulentViscosity','dragTerm','turbStress']
	if twoSize:
		globalParam =['qsMean','phiPart','vxPart','vyPart','vzPart','vxFluid','turbulentViscosity','dragTerm','turbStress', 'phiPart1','vxPart1','vyPart1','vzPart1','phiPart2','vxPart2','vyPart2','vzPart2']
	Save(path, globalParam)
	fileNumber +=1
	
#PHF	if twoSize!=1 and O.time>=150:
#	if O.time>=150:
#		print '\n End of the simulation, simulated 150s as required !\n '
#		O.pause()
#		O.exitNoBacktrace()

	#For reload, save sometimes the variables necessary to reload. (does not work anymore...)
	if (fileNumber*measurePeriod)%savePeriod == 0:
		#Save the simulation of yade the global python parameter in order to be able to recover the simulation
		O.save(scriptPath + '/'+'save.xml')
		# Save the global parameters which shouldn't be initialized to zero at the reloading. 
		globalParam = ['phiPart','fileNumber','vxFluid','turbulentViscosity','waterDepth','turbStress']
		Save(scriptPath + '/'+'paramSave.py', globalParam)


	##
	# Evaluation and printing of the Shields number or qsMean.

	if imposedFluidPro==1 or initImposedFluidPro==1:
		print 'qsMean', qsMean
	else:
		if twoSize:
			#Consider the average diameter in the layer in motion (which is the one to consider for Shields number)
			dummy = np.where(phiPart[:ndimz]<0.8*phiPartMax)
			phi1Sum = sum(phiPart1[dummy])
			phi2Sum = sum(phiPart2[dummy])
			#Volume weighted average diameter: N1Vp1 d1 + N2Vp2 d2/(N1Vp1 + N2Vp2) (*VCell/VCell)
			diameterMix = (diameterPart1*phi1Sum + diameterPart2*phi2Sum)/(phi1Sum + phi2Sum)
		else:
			#Only one size in the simulation...
			diameterMix = diameterPart
		shieldsNumber = densFluid*max(turbStress)/((densPart-densFluid)*diameterMix*abs(gravityVector[2]))
		print 'Shields number', shieldsNumber



###############
#########################################


##########		 		 		 		 		 	      ##########
### EVALUATE THE MEAN VOLUME FRACTION, DRAG AND STREAMWISE VELOCITY OF THE PARTICLES IN FUNCTION OF THE DEPTH###
##########		 		 		 		 		 	      ##########
#Arguments : ndimz, dz, length, Vlayer, groundPosition, radius, phiPart
#Modify : phiPart, vxPart, averageDrag
def averageProfile():
	global phiPart
	global averageDrag
	global vxPart
	global vyPart
	global vzPart
	global qsMean

	global phiPart1
	global averageDrag1
	global phiPart2
	global averageDrag2

	#Initialization of the concentration and the streamwise velocity. (two times the water height in order to take into account the particles which are also above the fluid in the concentration and velocity profile
	phiPart_loc = np.zeros(ndimz)		#Mean particle concentration in function of depth
	vxPart_loc = np.zeros(ndimz)		#Mean particle velocity in function of depth
	vyPart_loc = np.zeros(ndimz)		#Mean particle velocity in function of depth
	vzPart_loc = np.zeros(ndimz)		#Mean particle velocity in function of depth
	averageDrag_loc = np.zeros(ndimz)	#Mean drag

	phiPart1_loc = np.zeros(ndimz)
	averageDrag1_loc = np.zeros(ndimz)
	phiPart2_loc = np.zeros(ndimz)
	averageDrag2_loc = np.zeros(ndimz)
	if debugStress==1:
		global vxPartBIS
		global phiPartBIS
		global averageDragBIS
		global dragTermBIS
		averageDragBIS_loc = np.zeros(ndimz)
		phiPartBIS_loc = np.zeros(ndimz)
		vxPartBIS_loc = np.zeros(ndimz)
		dragTermBIS_loc = np.zeros(ndimz)


	#Import the vector of velocity fluctuation associated to each particle in C++ code
	VxFluct = addFluidForce.vFluctX
	VzFluct = addFluidForce.vFluctZ

        shiftPosition = 1.2*diameterPart        #To shift the position
in order to avoid to have negative position (not compatible with the
volume calculation used later)

	## Loop over the file to determine the concentration and the mean velocity
	for b in O.bodies:
		if isinstance(b.shape,Sphere):   # If it is a sphere
			z = b.state.pos[2]-groundPosition 	# Altitude of the particle redefined with the position of the bottom wall as the 
			radius = b.shape.radius
			Np = int(z/dz)			#Layer corresponding to the position of the center of the particle

			#Evaluate the drag force undergone by the particle (Dallavalle formulation + Richardson Zaki correction)
			uRel = Vector3(0.0,.0,0.0)
			if Np<ndimz and Np>=0 : uRel =  Vector3(vxFluid[Np]+VxFluct[b.id],0.0,VzFluct[b.id]) - b.state.vel
			fDrag = 0.5*pi*pow(radius,2)*densFluid*(0.44*uRel.norm()+24.4*kinematicViscoFluid/(2*radius))*pow(1-phiPart[Np],-expoRZ)*uRel

			if debugStress==1:
				dragTermBIS_loc[Np]+=fDrag[0]
				phiPartBIS_loc[Np]+=1.0
				vxPartBIS_loc[Np]+=b.state.vel[0]

			#Shift the position to avoid z<0 which induces some problems in the volume calculation. This is done here in order not to change the value of the drag. Evaluate all the parameters after this shift. (Np in particular is reevaluated
			z+=shiftPosition	
			Np = int(z/dz)			#new evaluation of Np after the shift. DO NOT REMOVE, IMPORTANT
			minZ = int((z-radius)/dz)	#Lowest layer occupied by the particle
			maxZ = int((z+radius)/dz)	#Highest layer occupied by the particle
			deltaCenter = z - Np*dz		#length between the real center of the particle and the position of the layer in which it is considered to be (Np).

			#loop from the bottom to the top layer inside which the particle is contained
			numLayer = minZ
			while numLayer <= maxZ:
				if numLayer>=0 and numLayer<ndimz:
					zInf = (numLayer-Np-1)*dz + deltaCenter
					zSup = (numLayer-Np)*dz + deltaCenter

					if numLayer==minZ: zInf = -radius
					if numLayer==maxZ: zSup = radius

					#evaluate the volume of the part of the particle contained in the layer considered (the formula results from the analytical resolution of the integral of a layer of sphere in cylindrical coordinate)
					V = math.pi*radius*radius*(zSup-zInf +(pow(zInf,3)-pow(zSup,3))/(3*radius*radius))

					#Add the quantities to the (not normalized for now) mean particle velocity and concentration
					phiPart_loc[numLayer]+=V
					vxPart_loc[numLayer]+=V*b.state.vel[0]
					vyPart_loc[numLayer]+=V*b.state.vel[1]
					vzPart_loc[numLayer]+=V*b.state.vel[2]
					averageDrag_loc[numLayer]+=V*fDrag[0]
					if twoSize:
						if radius==diameterPart1/2.0:
							phiPart1_loc[numLayer]+=V
							averageDrag1_loc[numLayer]+=V*fDrag[0]
						elif radius==diameterPart2/2.0:
							phiPart2_loc[numLayer]+=V
							averageDrag2_loc[numLayer]+=V*fDrag[0]
			 	#Increment for the next loop
				numLayer+=1

	#Normalize the weighted velocity/different mean quantities of each layer by the total particle volume in this layer
	n = 0
	for i in phiPart_loc:
		if i:  #if non zero, to avoid division by zero
			vxPart_loc[n]/=i
			vyPart_loc[n]/=i
			vzPart_loc[n]/=i
			averageDrag_loc[n]/=i
			if twoSize:
				if phiPart1_loc[n]:
					averageDrag1_loc[n]/=phiPart1_loc[n]
				if phiPart2_loc[n]:
					averageDrag2_loc[n]/=phiPart2_loc[n]

		n += 1
	#Normalization of the concentration
	phiPart_loc/=VLayer
	phiPart1_loc/=VLayer
	phiPart2_loc/=VLayer
	
	if debugStress==1:
		#Normalization
		n = 0
		for i in phiPartBIS_loc:#Loop over the number of part contained in each layer
			#Normalize the average drag and particle velocity
			if i:
				averageDragBIS_loc[n] = dragTermBIS_loc[n]/i
				vxPartBIS_loc[n]/=i
			n+=1
		phiPartBIS_loc*=(partVolume/VLayer)
		dragTermBIS_loc/=VLayer
		#Report the value in global variables
		vxPartBIS = vxPartBIS_loc
		phiPartBIS = phiPartBIS_loc
		averageDragBIS = averageDragBIS_loc 
		dragTermBIS = dragTermBIS_loc


	#Resize the vector to remove the effect of the shift previously introduced
	phiPart_loc = Resize(phiPart_loc,int(shiftPosition/dz))
	vxPart_loc = Resize(vxPart_loc,int(shiftPosition/dz))
	vyPart_loc = Resize(vyPart_loc,int(shiftPosition/dz))
	vzPart_loc = Resize(vzPart_loc,int(shiftPosition/dz))
	averageDrag_loc = Resize(averageDrag_loc,int(shiftPosition/dz))

	phiPart1_loc = Resize(phiPart1_loc,int(shiftPosition/dz))
	averageDrag1_loc = Resize(averageDrag1_loc,int(shiftPosition/dz))
	phiPart2_loc = Resize(phiPart2_loc,int(shiftPosition/dz))
	averageDrag2_loc = Resize(averageDrag2_loc,int(shiftPosition/dz))
	#Impose variables at the bottom to be consistent with the fluid formulation 
	averageDrag_loc[0] = 0.0
	vxPart_loc[0] = 0.0

	###Actualize the global variables
	phiPart = phiPart_loc
	vxPart =  vxPart_loc
	vyPart =  vyPart_loc
	vzPart =  vzPart_loc
	averageDrag = averageDrag_loc

	phiPart1 = phiPart1_loc
	averageDrag1 =  averageDrag1_loc
	phiPart2 = phiPart2_loc
	averageDrag2 =  averageDrag2_loc
	#And the one associated to the fluid force application
	addFluidForce.phiPart = phiPart_loc

	return
###############
#########################################



######	      				              ######
## Function to resize the vector removing the n first case##
######					              ######
def Resize(a,n):
	import numpy as np
	l = len(a)
	b = np.zeros(l)
	k = 0
	for i in range(n,l):
		b[k] = a[i]
		k+=1
	return b 
###############
#########################################


######					  ######
###  	    Save global variables	     ###
######					  ######

#Function to save global variables in a python file which can be re-executed for later
def Save(filePathName, globalVariables):
	f = open(filePathName,'w')
	f.write('from numpy import *\n')
	for i in globalVariables:
		f.write(i + ' = '+repr(globals()[i]) + '\n')
	f.close()
###############
#########################################
	

######					  ######
###    python Stress field calculation 	     ###
######					  ######
def stressFieldCalculation(ndimz,vCell,vxPart,dz):
	#Initialization
	granularTemp_loc = np.zeros(ndimz)
	N_loc =  np.zeros(ndimz)
	stressTensorField_loc = []
	kineticStressTensorField_loc = []
	for i in range(0,ndimz):
		stressTensorField_loc += [Matrix3(0,0,0,0,0,0,0,0,0)]
		kineticStressTensorField_loc += [Matrix3(0,0,0,0,0,0,0,0,0)]
	###
	## Dynamical stress tensor part
	###
	for b in O.bodies:
		Np = int((b.state.pos[2] - groundPosition)/dz)
		if Np>=0 and Np<ndimz:
			vFluct = b.state.vel - Vector3(vxPart[Np],vyPart[Np],vzPart[Np])
			stressTensorField_loc[Np]+= -1/vCell*b.state.mass*vFluct.outer(vFluct)
			kineticStressTensorField_loc[Np]+= -1/vCell*b.state.mass*vFluct.outer(vFluct)
			granularTemp_loc[Np] += 1/3.0*(pow(vFluct[0],2) + pow(vFluct[1],2) + pow(vFluct[2],2))
			N_loc[Np]+=1.0
	###
	## Love-Weber stress tensor part (contact
	###
	for Int in O.interactions:
		if O.bodies[Int.id1].dynamic or O.bodies[Int.id2].dynamic:
			Np1 = int((O.bodies[Int.id1].state.pos[2] - groundPosition)/dz)
			Np2 = int((O.bodies[Int.id2].state.pos[2] - groundPosition)/dz)
			#Vector joining the two particles (from id2 to id1)
			x12 = O.bodies[Int.id1].state.pos - O.bodies[Int.id2].state.pos
			#Contact vector of the two particles (defined as applied by id1 on id2)
			fContact = Int.phys.normalForce + Int.phys.shearForce
			#Remove the artificial effect of the periodic BC if one particle goes to the other side
			x12 -= O.cell.hSize*Int.cellDist
			#No reason to have one upon another, so check and define the min and max correpondingly
			if Np2==Np1:	#If the contact vector is entirely inside the layer, add it directly. 
				if Np1>=0 and Np1<ndimz:
					stressTensorField_loc[Np1]+= 1/vCell*fContact.outer(x12)	#outer product of contact force and branch vector		
			else:
				if Np1>Np2:	#Particle 1 above 2 (z component)
					minZ = Np2
					minPosZ = O.bodies[Int.id2].state.pos[2] - groundPosition
					maxZ = Np1
					maxPosZ = O.bodies[Int.id1].state.pos[2] - groundPosition
				elif Np2>Np1:	#Particle 2 above 1 (z component)
					minZ = Np1
					minPosZ = O.bodies[Int.id1].state.pos[2] - groundPosition
					maxZ = Np2
					maxPosZ = O.bodies[Int.id2].state.pos[2] - groundPosition

				#Loop over the layers containing the vector connecting the two particles
				numLayer = minZ
				while numLayer<=maxZ:
					if numLayer>=0 and numLayer<ndimz: 
						deltaZ = dz
						if numLayer==minZ:
                                                       deltaZ = dz - (minPosZ - minZ*dz)
						elif numLayer==maxZ:
                                                       deltaZ = maxPosZ - maxZ*dz

						branchVectCell=deltaZ/abs(x12[2])*Vector3(x12[0],x12[1],x12[2])
						#Add to stress tensor of the considered layer the part considered
                                                stressTensorField_loc[numLayer]+= 1/vCell*fContact.outer(branchVectCell)    #outer product
					#Increment the layer
					numLayer+=1
	#Normalization
	for n in range(0,len(N_loc)):
		if N_loc[n]>0:
			granularTemp_loc[n]/=N_loc[n]

        return
[stressTensorField_loc,kineticStressTensorField_loc,granularTemp_loc]


###############
#########################################


###DEBUG/COMPARISON WITH getStress
#To obtain exactly the same as getStress, you need to shift  groundPosition -> groundPosition - diameterPart + to remove the condition on the dynamics of the particles, to avoid the non contribution of the bottom wall and particles.
#This function might be also useful, it is reproducing exactly getStress in python
def getStressBIS(vCell):
	stressTensorField = Matrix3(0,0,0,0,0,0,0,0,0)
	for Int in O.interactions:
		x12 = O.bodies[Int.id1].state.pos - O.bodies[Int.id2].state.pos
		x12 -= O.cell.hSize*Int.cellDist
		fContact = Int.phys.normalForce + Int.phys.shearForce
		stressTensorField+= fContact.outer(x12)	#outer product of contact force a
	return stressTensorField/vCell
###############
#########################################

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.