//+
Solver.AutoCheck = 0;
SetFactory("OpenCASCADE");
//Mesh.CharacteristicLengthMin = 0.05;
//Mesh.CharacteristicLengthMax = 0.06;
Mesh.MshFileVersion = 2.2;
Mesh.CharacteristicLengthMin = 0.5;
Mesh.CharacteristicLengthMax = 0.5;

// Parameters
d_epi = 3.77;
d_endo = 3.75;
l_epi = 0.69;
l_endo = 0.49;
lr_epi = 1.1;
lr_endo = 1.05;
mu_base = 120.0/180.8*Pi;

a_epi = d_epi*Cosh(l_epi);
b_epi = d_epi*Sinh(l_epi);
c_epi = d_epi*Sinh(lr_epi);
a_endo = d_endo*Cosh(l_endo);
b_endo = d_endo*Sinh(l_endo);
c_endo = d_epi*Sinh(lr_endo);
ra_endo = 0.5*(a_epi+a_endo);
rb_endo = 0.5*(b_epi+b_endo);

// create entire left ventricle sphere
lvs_epi = newv;
Sphere(lvs_epi) = {0.0, 0.0, 0.0, 1.0, Pi/2-mu_base, Pi/2, 2*Pi};
// scale to ellipsoid reprenting left ventricle
lv_bl_tiss = newv;
lv_bl_tiss = Dilate{{0.,0.,0.},{b_epi,b_epi,a_epi}}{Volume{lvs_epi};};
////+
// create left ventricle blood sphere
lvs_endo = newv;
Sphere(lvs_endo) = {0.0, 0.0, 0.0, 1.0, Pi/2-1.1*mu_base, Pi/2, 2*Pi};
// scale to ellipsoid reprenting left ventriclar blood mass
lv_bl = newv;
lv_bl = Dilate{{0.,0.,0.},{b_endo,b_endo,a_endo}}{Volume{lvs_endo};};
//+
// create right ventricle half sphere
rvs_epi = newv;
Sphere(rvs_epi) = {0.0, 0.0, 0.0, 1.0, Pi/2-mu_base, Pi/2, Pi};
// scale to ellipsoid reprenting right ventricle
rv_bl_tiss = newv;
rv_bl_tiss = Dilate{{0.,0.,0.},{b_epi,c_epi,a_epi}}{Volume{rvs_epi};};
//+
// create right ventriclar blood half sphere
rvs_endo = newv;
Sphere(rvs_endo) = {0.0, 0.0, 0.0, 1.0, Pi/2-1.1*mu_base, Pi/2, Pi};
// scale to ellipsoid reprenting right ventriclar blood mass
rv_bl = newv;
rv_bl = Dilate{{0.,0.,0.},{rb_endo,c_endo,ra_endo}}{Volume{rvs_endo};};

// create lv tissue region
lv_tiss = newv;
BooleanDifference(lv_tiss) = {Volume{lv_bl_tiss}; Delete;}{Volume{lv_bl}; Delete;};

// create rv tissue region
rv_tiss = newv;
BooleanDifference(rv_tiss) = {Volume{rv_bl_tiss}; Delete;}{Volume{rv_bl}; Delete;};

// Determine rv tissue region different from lv tissue region
rv_tiss_fin = newv;
BooleanDifference(rv_tiss_fin) = {Volume{rv_tiss}; Delete;}{Volume{lv_tiss};};
// merge lv and rv tissue regions to create tissue region for heart
tissue = newv;
BooleanUnion(tissue) = {Volume{lv_tiss}; Delete; }{ Volume{rv_tiss_fin}; Delete; };
