# 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
# Transformation of a square of ethanol into a circle due to the tension surface force exerted by the air surrounding the square

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)

#physical parameters
g = np.array([0,0])               
rho0 = 1.117                                            # air density
rho1 = 785.92                                           # ethanol density
nu0 = 1.57e-5                                           # air kinematic viscosity
nu1 = 1.2e-3/rho1                                       # ethanol kinematic viscosity
R = 0.01*np.pi**0.5/2*10                                # square lenght
sigma = 22.39e-3                                        # surface tension

#numerical parameters
tEnd = 200                                              # final time
dt = 1e-3                                               # time step
outf = 10                                               # number of iterations between output files
ii = 0                                                  # initial iteration
t = 0                                                   # initial time

#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1],sigma=sigma)
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
fluid.set_wall_boundary("Top",pressure=0)
fluid.set_strong_boundary("Top",0,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Bottom",0,0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Lateral",0,0)
fluid.set_strong_boundary("Lateral",1,0)

# Solution fields and coordinates
s = fluid.solution()
c = np.ndarray((fluid.n_nodes()))
x = fluid.coordinates()

# Initial concentration field
c[:] = 0
c[np.logical_and(np.abs(x[:,0])<R,np.abs(x[:,1])<R)] = 1
fluid.set_concentration_cg(c)
fluid.write_vtk(outputdir,0,0)


tic = time.time()
#
# COMPUTATION LOOP
#
while t < tEnd : 
    time_integration.iterate(fluid,None,dt)
    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