Include "PlaneWaveBiperiodic3D_data.geo";
SetFactory("OpenCASCADE");
lc_super = lc/Sqrt[epsr_super_re];
lc_subs  = lc/Sqrt[epsr_subs_re];

// domains, bottom to top
Box(1) = {-d/2, -d/2, -h_subs-h_pml_subs, d , d , h_pml_subs};
Box(2) = {-d/2, -d/2, -h_subs, d , d , h_subs};
Box(3) = {-d/2, -d/2, 0, d , d , h_super};
Box(4) = {-d/2, -d/2, h_super, d , d , h_pml_super};
Coherence;

// get Bloch periodic faces
e=1*nm;
BlochXp() = Surface In BoundingBox{  d/2-e, -d/2-e, -h_subs-h_pml_subs-e,  d/2+e,  d/2+e, h_super+h_pml_super+e};
BlochXm() = Surface In BoundingBox{ -d/2-e, -d/2-e, -h_subs-h_pml_subs-e, -d/2+e,  d/2+e, h_super+h_pml_super+e};
BlochYm() = Surface In BoundingBox{ -d/2-e, -d/2-e, -h_subs-h_pml_subs-e,  d/2+e, -d/2+e, h_super+h_pml_super+e};
BlochYp() = Surface In BoundingBox{ -d/2-e,  d/2-e, -h_subs-h_pml_subs-e,  d/2+e,  d/2+e, h_super+h_pml_super+e};
Periodic Surface{BlochXp()} = {BlochXm()} Translate{ d , 0 , 0 };
Periodic Surface{BlochYp()} = {BlochYm()} Translate{ 0 , d , 0 };

// get top / bottom PML/media interfaces
Sbot() = Surface In BoundingBox{-d/2-e, -d/2-e, -h_subs-e , d/2+e, d/2+e, -h_subs+e };
Stop() = Surface In BoundingBox{-d/2-e, -d/2-e,  h_super-e, d/2+e, d/2+e,  h_super+e};


Physical Volume("PMLBOT",1) = {1};
Physical Volume("SUBS"  ,2) = {2};
Physical Volume("SUPER" ,3) = {3};
Physical Volume("PMLTOP",4) = {4};
Physical Surface("BXM" ,101 ) = BlochXm();
Physical Surface("BXP" ,102 ) = BlochXp();
Physical Surface("BYM" ,201 ) = BlochYm();
Physical Surface("BYP" ,202 ) = BlochYp();
Physical Surface("STOP" ,301 ) = Stop();
Physical Surface("SBOT" ,302 ) = Sbot();

Characteristic Length{PointsOf{Physical Volume{4};}} = lc_super;
Characteristic Length{PointsOf{Physical Volume{3};}} = lc_super;
Characteristic Length{PointsOf{Physical Volume{2};}} = lc_subs;
Characteristic Length{PointsOf{Physical Volume{1};}} = lc_subs;
