//+
SetFactory("OpenCASCADE");
Mesh.MshFileVersion = 2.2;
Mesh.CharacteristicLengthMin = 0.2;
Mesh.CharacteristicLengthMax = 0.2;
sph1 = newv;
Sphere(sph1) = {0.0, 0.0, 0.0, 1.0, 0.0, Pi/2, 2*Pi};
sph2 = newv;
Sphere(sph2) = {0.0, 0.0, 0.0, 2.0, 0.0, Pi/2, 2*Pi};

// theta range from th_l to th_u
th_ld = 10.0;
th_ud = 40.0;

// offset in phi direction
phioffd = 160.0;

// phi range from phioff-phid to phioff+phid
phid = 20.0;

// theta range from th_l to th_u
th_l = th_ld*Pi/180.0;
th_u = th_ud*Pi/180.0;
th_c = 0.5*(th_l+th_u);

// offset in phi direction
phioff = phioffd*Pi/180.0;

// phi range from phioff-phi to phioff+phi
phi = phid*Pi/180.0;


origin = newp;
Point(origin) = {0.0,0.0,0.0};
shell = newv;
BooleanDifference(shell) = {Volume{sph2};}{Volume{sph1}; };

//// parameter for high mesh density near sub region
md = 0.05;

ll = newp;
Point(ll) = {Cos(th_l)*Cos(phioff-phi),Cos(th_l)*Sin(phioff-phi),Sin(th_l),md};
ul = newp;
Point(ul) = {Cos(th_u)*Cos(phioff-phi),Cos(th_u)*Sin(phioff-phi),Sin(th_u),md};
ur = newp;
Point(ur) = {Cos(th_u)*Cos(phioff+phi),Cos(th_u)*Sin(phioff+phi),Sin(th_u),md};
lr = newp;
Point(lr) = {Cos(th_l)*Cos(phioff+phi),Cos(th_l)*Sin(phioff+phi),Sin(th_l),md};

centrept = newp;
Point(centrept) = {Cos(th_c)*Cos(phioff),Cos(th_c)*Sin(phioff),Sin(th_c),md};
//
//Point{ll,ul,ur,lr,centrept} In Volume {shell};
//// create circular arcs for patch sides
c1 = newl;
Circle(c1) = {ll,origin,ul};
c2 = newl;
Circle(c2) = {ul,origin,ur};
c3 = newl;
Circle(c3) = {ur,origin,lr};
c4 = newl;
Circle(c4) = {lr,origin,ll};
//Curve{c1,c2,c4,c4} In Volume{shell};
//
//// create curve loop and surface
cl1 = newll;
Curve Loop(cl1) = {c1,c2,c3,c4};
patch = news;
Surface(patch) = {cl1};

//Surface{patch} In Volume{shell};
//
//// get coordinates of centre point of patch
cpoint[] = Point{centrept};
//
//// define normal vector to surface at centre point
norm_x = 2.0*cpoint[0];
norm_y = 2.0*cpoint[1];
norm_z = 2.0*cpoint[2];
//
//// get length of vector
len = Sqrt(norm_x^2+norm_y^2+norm_z^2);
//
//// define unit vector scaled by thickness of ischaemic region.
thickness = 0.5;
u_norm_x = thickness*norm_x/len;
u_norm_y = thickness*norm_y/len;
u_norm_z = thickness*norm_z/len;
//
//// create sub region volume
sub_vol() = Extrude {u_norm_x,u_norm_y,u_norm_z} {Surface{patch}; };
//
//////////// define two regions
region1() = BooleanDifference {Volume{shell};}{Volume{sub_vol[1]};};
region2() = BooleanIntersection {Volume{shell};Delete;}{Volume{sub_vol[1]}; Delete;};
//
Printf('region 1: ',region1());
Printf('region 2: ',region1());
//
//// Idea from Fossies t16.geo (tutorial 16)
r1_surf() = Unique(Abs(Boundary{ Volume{region1}; }));
r2_surf() = Unique(Abs(Boundary{ Volume{region2}; }));
Coherence;
//
//Mesh 2;
//
Printf('Region 1 Surfaces:  ',r1_surf());
Printf('Region 2 Surfaces:  ',r2_surf());
