# MigFlow - Copyright (C) <2010-2020>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
# 	
# This program (MigFlow) is free software: 
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
# 
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files).  If not, 
# see <http://www.gnu.org/licenses/>.

#!/usr/bin/env python

# TESTCASE DESCRIPTION
# Study of the drag on a fixed cylinder in a 3D immersed granular flow
from migflow import scontact
from migflow import fluid
from migflow import time_integration
import numpy as np
import os
import time
import shutil
import random
import sys
listeVit = [14,36,64,95,130,169,211,255,303,354,408,467]
inter = int(sys.argv[1])
index = int(inter/5)
numDep = inter%5
# Physical parameters
vit = -0.001*listeVit[index]                   #stream velocity
muwall = 0.3                                   #friction coefficient between the walls and the particles
mupart = 0.15                                  #friction coefficient between particles
rho = 1000                                     #fluid density
nu = 1e-6                                      #fluid kinematic viscosity
g = np.array([0,-9.81,0])                      #gravity
# Numerical parameters
dt = 1e-3                                      #time step
t = 0                                          #initial time
tEnd = 30                                      #final time
ii = 0                                         #iteration number
outf = 2500                                    #iterations between data frames  
# Define output directory
outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)
filename = outputdir + "/Drag.csv"
# 
# PARTICLE PROBLEM
#
#Injected particles
p = scontact.ParticleProblem(3,True,True)
p.read_vtk("depot/"+ str(numDep),2)
p.remove_particles_flag((p.position()[:,1]-p.r()[:,0]>0.5)*(p.position()[:,1]+p.r()[:,0] <0.52))
p.velocity()[:,1] = vit
#Flow particles
p2 = scontact.ParticleProblem(3,True,True)
p.read_vtk("depot/"+ str(numDep),2)
p2.remove_particles_flag(p2.position()[:,1]+p2.r()[:,0] <0.5)
p2.clear_boundaries()
p2.load_msh_boundaries("mesh.msh", ["Inner","Left", "Right","Front","Rear"],material = "Steel")
p2.set_friction_coefficient(mupart,"Sand","Sand")
p2.set_friction_coefficient(muwall,"Sand","Steel")
p2.write_vtk(outputdir, 0, 0)
#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(3,g,[nu*rho],[rho],drag_in_stab=0)
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Inner")
fluid.set_wall_boundary("Left", velocity = [0,0,0])
fluid.set_wall_boundary("Right", velocity = [0,0,0])
fluid.set_open_boundary("Bottom", velocity = [0,vit,0])
fluid.set_wall_boundary("Front", velocity = [0,0,0])
fluid.set_wall_boundary("Rear", velocity = [0,0,0])
fluid.set_open_boundary("Top",pressure = 0)
fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces())
fluid.export_vtk(outputdir, 0,0)
# 
# COMPUTATION LOOP
#
def accumulate(F,n_divide):
  F+=np.sum(p2.get_boundary_forces("Inner"),axis=0)/n_divide
while t < tEnd : 
     p2.forced_flag()[(p2.position()[:,1] + p2.r()[:,0] <= -0.5)] = 1
     p2.velocity()[p2.forced_flag()==1,:] = [0, vit, 0]
     p2.omega()[p2.forced_flag()==1,:] = [0,0,0]
     if  max(p2.position()[:,1]+p2.r()[:,0]) < 0.5 and t < 21 :
       p2.add_particles(p.position(),p.r(),p.mass(),v=p.velocity(),tag="Sand",forced=p.forced_flag(),contact_forces=p.contact_forces())
     p2.remove_particles_flag( (p2.position()[:,1] + p2.r()[:,0] >-0.52))
     F = np.zeros(3)
     fluid.set_particles(p2.mass(), p2.volume(), p2.position(), p2.velocity(), p2.contact_forces())
     time_integration.iterate(None,p2,dt,1,contact_tol=1e-6,external_particles_forces=g*p2.mass(),after_sub_iter=lambda n_divide : accumulate(F,n_divide))
     with open(filename,"a") as file1:
        F += fluid.boundary_force()[0,:]
        file1.write( str(F[0])+";"+str(F[1])+";"+str(F[2])+";"+str(t)+"\n")
     t += dt   
     # Output files writing
     if ii %outf == 0 :
         ioutput = int(ii/outf) + 1
         p2.write_vtk(outputdir, ioutput, t)
         fluid.write_vtk(outputdir, ioutput, t)       
     ii += 1
     print("%i : %.2g/%.2g" % (ii, t, tEnd))