import gmsh
import sys
import numpy as np

# set up parameters (eventually will read from file). 
# Data from streeter (1977?)
d_epi = 3.77
d_endo = 3.75
l_epi = 0.69
l_endo = 0.49
mu_base = 120.0/180.8*np.pi

# create unit spheres
radius = 1.0

# set CharacteristicLengthMax
char_l_max = 0.2
#char_l_max = 0.06

# my parameters to create the RV
lr_epi = 1.1
lr_endo = 1.05

# derived quantities
a_epi = d_epi*np.cosh(l_epi)
b_epi = d_epi*np.sinh(l_epi)
c_epi = d_epi*np.sinh(lr_epi)
a_endo = d_endo*np.cosh(l_endo)
b_endo = d_endo*np.sinh(l_endo)
c_endo = d_epi*np.sinh(lr_endo)
ra_endo = 0.5*(a_epi+a_endo)
rb_endo = 0.5*(b_epi+b_endo)

# start of main code
gmsh.initialize(sys.argv)
gmsh.model.add("streeter")

# to print errors
gmsh.option.setNumber("General.Terminal", 1)
gmsh.option.setNumber("Geometry.ToleranceBoolean",1.e-4)
gmsh.option.setNumber("Mesh.MshFileVersion",2.2)

# mesh parameters
gmsh.option.setNumber('Mesh.CharacteristicLengthMax',char_l_max)

# create spheres for left and right blood and blood+tissue regions
# create entire left ventricle
lv_bl = gmsh.model.occ.addSphere(0,0,0, radius, -1, np.pi/2-mu_base,
                                   np.pi/2, 2*np.pi)
gmsh.model.occ.dilate([(3,lv_bl)],0.,0.,0.,b_epi,b_epi,a_epi)

#create left ventricle blood region
lv_bl_tiss = gmsh.model.occ.addSphere(0,0,0, radius, -1, np.pi/2-1.1*mu_base,
                                   np.pi/2, 2*np.pi)
gmsh.model.occ.dilate([(3,lv_bl_tiss)],0.,0.,0.,b_endo,b_endo,a_endo)

# create right ventricle (half ellipsoid)
rv_bl = gmsh.model.occ.addSphere(0,0,0, radius,-1, np.pi/2-mu_base,
                                   np.pi/2, np.pi)
gmsh.model.occ.dilate([(3,rv_bl)],0.,0.,0.,b_epi,c_epi,a_epi)

# create right ventricle blood region
rv_bl_tiss = gmsh.model.occ.addSphere(0,0,0, radius,-1, np.pi/2-1.1*mu_base,
                                   np.pi/2, np.pi)
gmsh.model.occ.dilate([(3,rv_bl_tiss)],0.,0.,0.,rb_endo,c_endo,ra_endo)

# get lv tissue region
lv_tiss, _ = gmsh.model.occ.cut([(3,lv_bl)],[(3,lv_bl_tiss)], -1, True, True)
print('lv_tiss = ',lv_tiss)
# get rv tissue region
rv_tiss, _  = gmsh.model.occ.cut([(3,rv_bl)],[(3,rv_bl_tiss)], -1, True, True)
print('rv_tiss = ',rv_tiss)
# get rv tissue region not in lv
rv_tiss_fin, _  = gmsh.model.occ.cut([rv_tiss],[lv_tiss], -1, True, False)
print('rv_tiss_fin = ',rv_tiss_fin)
gmsh.model.occ.synchronize()
# merge lv and rv tissue regions 
tissue, _  = gmsh.model.occ.fuse([lv_tiss],[rv_tiss_fin], -1, True, True)
gmsh.model.occ.synchronize()
#
## uncomment the following line to display the mesh
gmsh.fltk.run()
gmsh.finalize()