import espressomd
import espressomd.lb

from espressomd import lbboundaries
from espressomd import shapes
from espressomd import interactions

import object_in_fluid as oif
import numpy as np
import os, sys
import random

simNo = sys.argv[1]
os.mkdir("output/sim" + str(simNo))

boxX = 40.0
boxY = 20.0
boxZ = 20.0
system = espressomd.System(box_l=[boxX, boxY, boxZ])
system.cell_system.skin = 0.2
time_step = 0.1
system.time_step = time_step

# creating the template for RBCs
type = oif.OifCellType(nodes_file="input/rbc374nodes.dat", triangles_file="input/rbc374triangles.dat",
                        check_orientation=False, system=system, ks=0.02, kb=0.02, kal=0.02,
                        kag=0.9, kv=0.5, normal=True, resize=[3.91, 3.91, 3.91])

# creating the RBC
cell = oif.OifCell(cell_type=type, particle_type=0, origin=[10.0,10.0,6.0], rotate=[np.pi/2,np.pi/2,0.0])

# fluid
lb_params = {'agrid': 1, 'dens': 1, 'visc': 1.5, 'tau': system.time_step, 'ext_force_density': [0.002, 0.0, 0.0]}
lbf = espressomd.lb.LBFluid(**lb_params)
system.actors.add(lbf)
system.thermostat.set_lb(LB_fluid=lbf, gamma=1.5)

maxCycle = 20
# main integration loop
cell.output_vtk_pos_folded(file_name="output/sim" + str(simNo) + "/cell_0.vtk")
integr_steps = 100
for i in range(1,maxCycle):

    flow_through_membrane = 0
    for t in cell.mesh.triangles:
        posA = t.A.get_pos()
        posB = t.B.get_pos()
        posC = t.C.get_pos()
        velDiffA = t.A.get_vel() - lbf.get_interpolated_velocity(posA)
        velDiffB = t.B.get_vel() - lbf.get_interpolated_velocity(posB)
        velDiffC = t.C.get_vel() - lbf.get_interpolated_velocity(posC)
        avg_vel = np.array(velDiffA + velDiffB + velDiffC)/3.0
        normal = np.array(oif.get_triangle_normal(posA, posB, posC))
        proj_avg_vel = (np.inner(avg_vel, normal)/np.inner(normal, normal)) * normal
        size_proj_avg_vel = oif.norm(proj_avg_vel)
        partial_volume = size_proj_avg_vel * oif.area_triangle(posA, posB, posC)
        flow_through_membrane = flow_through_membrane + partial_volume

    cell_volume = cell.volume()
    print "time: ", str(i * time_step * integr_steps), "flow: ", str(flow_through_membrane), "% volume: ", str(flow_through_membrane/cell_volume*100)
    system.integrator.run(steps=integr_steps)
    cell.output_vtk_pos_folded(file_name="output/sim" + str(simNo) + "/cell_" + str(i) + ".vtk")