# 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
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

outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

t = 0
ii = 0


#physical parameters
g =  0#-9.81                
rhof = 1030                                      # fluid density
rhop = 2450                                     # grains density
r=154e-6                                        # grains radii
R = 0.01                                      # drop radius
compacity = .20                                 # solid volume fraction in the drop
phi = 1 - compacity
nuf = 1.17/rhof                                   # kinematic viscosity
muf = nuf*rhof                                     # dynamic viscosity

##mixture properties
rhom = (1-phi)*rhop+phi*rhof                   #mixture density
num = nuf*(1+5/2*phi)                          #Einstein viscosity

print("rhom = %g, num = %g" % (rhom,num))


#numerical parameters
tEnd = 200                                     # final time
dt = 1e-2                                       # time step
outf = 5                                       # number of iterations between output files

def outerBndV(x) :
    print(.4*max(np.sin(t*np.pi*2./1),0))
    return 0.4*max(np.sin(t*np.pi*2./1),0)

fluid = fluid.FluidProblem(2,g,[nuf*rhof,num*rhom],[rhof,rhom],coeff_stab=1e-6)
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
fluid.set_wall_boundary("Top",pressure=0)
ii = 0
t = 0

s = fluid.solution()
c = np.ndarray((fluid.n_nodes()))
x = fluid.coordinates()

c[:] = 1
c[np.logical_and(np.abs(x[:,0])< R,np.abs(x[:,1])<R)] = 0
#for i in range(len(x[:,0])):
#    z = (x[i,0])**2+(x[i,1])**2
#    R1 = (0.7*R)**2
#    R2 = (1.3*R)**2
#    #s[i,1] = -0.01
#    c[i] = min(max(0,1/(R2-R1)*z-R1/(R2-R1)),1)

fluid.set_concentration_cg(c)


#set initial_condition

fluid.write_vtk(outputdir,0,0)

ii = 0
tic = time.time()
while t < tEnd : 
    #Fluid solver
    time_integration.iterate(fluid,None,dt)
    #if (ii%11==0 and ii != 0):
    #  fluid.adapt_mesh(1e-2,1e-4,5000)
    t += dt
    #Output files writting
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
        fluid.write_vtk(outputdir, ioutput, t)
    ii += 1