# 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 granular flow
from migflow import scontact
from migflow import fluid
from migflow import time_integration
import numpy as np
import os
def genInitialPosition(filename, r, rhop) :
    """Set all the particles centre positions and create the particles objects to add in the computing structure.
    
    Keyword arguments:
    filename -- name of the output file
    r -- mean radius of the particles
    ly -- domain height
    lx -- domain width 
    rhop -- particles density        
    """
    #Particles structure builder
    p = scontact.ParticleProblem(2, friction_enabled = True,rotation_enabled=True)
    #Load mesh.msh file specifying physical boundaries names
    p.load_msh_boundaries("mesh.msh", ["Inner"], material = "Steel")
    p.load_msh_boundaries("mesh.msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"], material = "Plexi")
    p.load_msh_boundaries("mesh.msh", ["LockDown"], material = "Steel")
    #Particles centre are placed on a regular grid
    # Add a grain at each centre position
    for y in np.arange(1+r, 2, 2*r*1.05):
            for x in np.arange(-0.295+r*1.05, 0.295-r*1.05, 2*r*1.05):
                r1 = np.random.uniform(.95*r, 1.05*r)
                p.add_particle((x, y), r1, r1**2 * np.pi * rhop,"Sand");  
    p.write_vtk(outputdir, 0, 0)
# Physical parameters
g = np.array([0,-9.81])
r = 1e-2
rhop = 2500
mupart= 0.12
# Numerical parameters
dt = 1e-3                                      #time step
t = 0                                          #initial time
tEnd = 13                                      #final time
ii = 0                                         #iteration number
outf = 100                                     #iterations between data frames  
# Define output directory
outputdir = "output" 
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)
#Creating file for drag data
filename = outputdir + "/Drag.csv"
# 
# PARTICLE PROBLEM
#
genInitialPosition(outputdir, r, rhop)
p = scontact.ParticleProblem(2,friction_enabled = True,rotation_enabled=True)
p2 = scontact.ParticleProblem(2,friction_enabled = True,rotation_enabled=True)
p.read_vtk(outputdir, 0)
p2.read_vtk(outputdir, 0)
x2 = np.copy(p2.position())
p2.remove_particles_flag((x2[:,1] + p2.r()[:,0] <1.15)*(x2[:,1] + p2.r()[:,0] >1.1)*(x2[:,0] + p2.r()[:,0] <0.15)*(x2[:,0] + p2.r()[:,0] >-0.15))
p2.position()[:,1] -= 0.6
p.clear_boundaries()
p.load_msh_boundaries("mesh.msh", ["Inner"], material = "Steel")
p.load_msh_boundaries("mesh.msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"], material = "Plexi")
p.load_msh_boundaries("mesh.msh", ["LockDown"], material = "Steel")
p.set_friction_coefficient(mupart,"Sand","Sand")
p.set_friction_coefficient(0.1,"Sand","Steel")
p.set_friction_coefficient(0.1,"Sand","Plexi")
p.write_vtk(outputdir, 0, 0)
# 
# COMPUTATION LOOP
#
opened = 0
def accumulate(F,ndivide):
  F+=np.sum(p.get_boundary_forces("Inner"),axis=0)/ndivide
while t < tEnd : 
    if t>=3 and opened == 0:
        opened = 1
        p.clear_boundaries()
        p.load_msh_boundaries("mesh.msh", ["Inner"],material="Steel")
        p.load_msh_boundaries("mesh.msh", ["Right","Left","RightUp","RightDown","LeftUp","LeftDown"],material="Plexi")
#Adding new particles
    if (max(p.position()[:,1])+max(p.r()) <0.5) and t > 3 and t < 25:
      p.add_particles(p2.position(),p2.r(),p2.mass(),v=p2.velocity(),tag="Sand",forced=p2.forced_flag(),contact_forces=p2.contact_forces())
#Solving the contacts
    F = np.zeros(2)
    time_integration.iterate(None,p,dt,min_nsub=1,contact_tol=1e-7,external_particles_forces=g*p.mass(),after_sub_iter= lambda ndivide : accumulate(F,ndivide))
#Removing particles outside the hopper
    p.remove_particles_flag( (p.position()[:,1] >-0.6))
#Computing mean velocity and writing to file
    with open(filename,"a") as file1:
        file1.write(str(F[0])+";"+str(F[1])+";"
                    +str(t)+"\n")
    t += dt
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        p.write_vtk(outputdir, ioutput, t)      
    ii += 1
    print("%i : %.2g/%.2g" % (ii, t, tEnd))