# 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
from migflow import lmgc90Interface
from pylmgc90 import pre
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(3, friction_enabled = True,rotation_enabled=True)
    #Load mesh.msh file specifying physical boundaries names
    p.load_msh_boundaries("mesh04" + ".msh", ["Inner","Right","Left","Front","Rear","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.7, 2*r*1.05):
            for x in np.arange(-0.19+r*1.05, 0.19-r*1.05, 2*r*1.05):
                for z in np.arange(0+r*1.05,0.05-1.05*r,2*r*1.05):
                    r1 = np.random.uniform(.95*r, 1.05*r)
                    p.add_particle((x, y,z), r1, (4/3)*r1**3 * np.pi * rhop,"Sand");  
    p.write_vtk(outputdir, 0, 0)
# Physical parameters
g = np.array([0,-9.81,0])
r = 1.1e-2
rhop = 2500
# Numerical parameters
dt = 1e-3                                      #time step
t = 0                                          #initial time
tEnd = 40                                      #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(3,friction_enabled = True,rotation_enabled=True)
p.set_friction_coefficient(0.3,"Sand","Sand")
p.set_friction_coefficient(0.3,"Sand","Steel")
p2 = scontact.ParticleProblem(3,friction_enabled = True,rotation_enabled=True)
p.read_vtk(outputdir, 0)
p2.read_vtk(outputdir, 0)
x2 = p2.position()
p2.remove_particles_flag((x2[:,1] + p2.r()[:,0] <1.25)*(x2[:,1] + p2.r()[:,0] >1.1)*(x2[:,0] + p2.r()[:,0] <0.18)*(x2[:,0] + p2.r()[:,0] >-0.18))
p2.position()[:,1] -= 0.5
p.clear_boundaries()
p.load_msh_boundaries("mesh" + ".msh", ["Inner","Right","Left","Front","Rear","LockDown"], material = "Steel")

p.write_vtk(outputdir, 0, 0)
# 
# COMPUTATION LOOP
#
opened = 0
def accumulate(F):
  F+=np.sum(p.get_boundary_forces("Inner"),axis=0)
while t < tEnd : 
    if t>=10 and opened == 0:
        opened = 1
        p.clear_boundaries()
        p.load_msh_boundaries("mesh" + ".msh", ["Inner","Right","Left","Front","Rear"],material="Steel")
#Adding new particles
    if (max(p.position()[:,1])+max(p.r()) <0.6) and t > 10 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(3)
    time_integration.iterate(None,p,dt,min_nsub=1,contact_tol=1e-6,external_particles_forces=g*p.mass(),after_sub_iter= lambda subdt : accumulate(F))
#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(F[2])+";"
                    +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))