// Sphere decomposed in concentric spheres
// For subdomain Omega_i, Every normal vector is outwardly directed EXCEPT
// for the normal on the scatterer which is inwardly directed
// Each layer has the same volume ; charachteristic length (lc) grows linearly
// from one interface to the other

Include "sphere_concentric_data.geo";

Solver.AutoMesh = -1; // the geometry script generates the mesh
Mesh.Optimize = 1;

If(N_DOM < 1)
  Printf("Hey Dude, stop kidding me :-)");
  Exit;
EndIf

Point(1) = {0,0,0,LC}; // center of the spheres

myLc = LC;

//average volume
vol_tot = 4/3*Pi*(R_EXT^3 - R_INT^3);
vol = vol_tot/N_DOM;

//Radius of sphere
R[0] = R_INT;

For cpt In {0:N_DOM}
  // r = (cpt+1)*R_INT; // constant subdomain radius -- problem size will depend on N_DOM

  r = R[cpt]; // constant subdomain volume

  R += (3./4./Pi*vol+R[cpt]^3.)^(1./3.); // compute next radius -- will be used at next step

  //scatterer (intern sphere)
  p_down = newp; 	Point(p_down) = {0, -r, 0, myLc};
  p_right = newp; Point(p_right) = {r, 0, 0, myLc};
  p_up = newp; 	Point(p_up) = {0, r, 0, myLc};
  p_left = newp; 	Point(p_left) = {-r, 0, 0, myLc};
  p_front = newp; Point(p_front) = {0, 0, r, myLc};
  p_back = newp; 	Point(p_back) = {0, 0, -r, myLc};

  l_dr = newl; Circle(l_dr) = {p_down, 1, p_right};
  l_ru = newl; Circle(l_ru) = {p_right, 1, p_up};
  l_ul = newl; Circle(l_ul) = {p_up, 1, p_left};
  l_ld = newl; Circle(l_ld) = {p_left, 1, p_down};

  l_uf = newl; Circle(l_uf) = {p_up, 1, p_front};
  l_fd = newl; Circle(l_fd) = {p_front, 1, p_down};
  l_db = newl; Circle(l_db) = {p_down, 1, p_back};
  l_bu = newl; Circle(l_bu) = {p_back, 1, p_up};

  l_lf = newl; Circle(l_lf) = {p_left, 1, p_front};
  l_fr = newl; Circle(l_fr) = {p_front, 1, p_right};
  l_rb = newl; Circle(l_rb) = {p_right, 1, p_back};
  l_bl = newl; Circle(l_bl) = {p_back, 1, p_left};

  ll_drf = newll; Line Loop(ll_drf) = {l_dr, -l_fr, l_fd};
  ll_urf = newll; Line Loop(ll_urf) = {l_ru, l_uf, l_fr};
  ll_ulf = newll; Line Loop(ll_ulf) = {l_ul, l_lf, -l_uf};
  ll_dlf = newll; Line Loop(ll_dlf) = {l_ld, -l_fd, -l_lf};

  ll_drb = newll; Line Loop(ll_drb) = {l_db, -l_rb, -l_dr};
  ll_urb = newll; Line Loop(ll_urb) = {l_bu, -l_ru, l_rb};
  ll_ulb = newll; Line Loop(ll_ulb) = {-l_bu, l_bl, -l_ul};
  ll_dlb = newll; Line Loop(ll_dlb) = {-l_db, -l_ld, -l_bl};

  surf_drf[cpt] = news; Ruled Surface(surf_drf[cpt]) = {ll_drf};
  surf_urf[cpt] = news; Ruled Surface(surf_urf[cpt]) = {ll_urf};
  surf_ulf[cpt] = news; Ruled Surface(surf_ulf[cpt]) = {ll_ulf};
  surf_dlf[cpt] = news; Ruled Surface(surf_dlf[cpt]) = {ll_dlf};

  surf_drb[cpt] = news; Ruled Surface(surf_drb[cpt]) = {ll_drb};
  surf_urb[cpt] = news; Ruled Surface(surf_urb[cpt]) = {ll_urb};
  surf_ulb[cpt] = news; Ruled Surface(surf_ulb[cpt]) = {ll_ulb};
  surf_dlb[cpt] = news; Ruled Surface(surf_dlb[cpt]) = {ll_dlb};

  surf_loop[cpt] = newsl; Surface Loop(surf_loop[cpt]) =
       {surf_drf[cpt], surf_urf[cpt], surf_ulf[cpt],
        surf_dlf[cpt], surf_drb[cpt], surf_urb[cpt],
        surf_ulb[cpt], surf_dlb[cpt]};

  If (cpt > 0)
    v = newv; vl[] += v; Volume(v) = {surf_loop[cpt-1], surf_loop[cpt]};
  EndIf
EndFor

If(StrCmp(OnelabAction, "check")) // only mesh if not in onelab check mode
  Mesh 3;
  CreateDir Str(DIR);

  For cpt In {1:N_DOM}
    Delete Physicals;

    If (cpt == 1)
      Physical Surface(-1001) = {surf_drf[cpt-1], surf_urf[cpt-1],
        surf_ulf[cpt-1], surf_dlf[cpt-1], surf_drb[cpt-1], surf_urb[cpt-1],
        surf_ulb[cpt-1], surf_dlb[cpt-1]};
      Physical Surface((4+cpt)*1000) = {surf_drf[cpt], surf_urf[cpt],
        surf_ulf[cpt], surf_dlf[cpt], surf_drb[cpt], surf_urb[cpt],
        surf_ulb[cpt], surf_dlb[cpt]};
      Physical Volume(6000+cpt) = {vl[cpt-1]};
    EndIf

    If (cpt == N_DOM)
      Physical Surface(-(3+cpt)*1000) = {surf_drf[cpt-1], surf_urf[cpt-1],
        surf_ulf[cpt-1], surf_dlf[cpt-1], surf_drb[cpt-1], surf_urb[cpt-1],
        surf_ulb[cpt-1], surf_dlb[cpt-1]};
      Physical Surface(2000+N_DOM) = {surf_drf[cpt], surf_urf[cpt], surf_ulf[cpt],
        surf_dlf[cpt], surf_drb[cpt], surf_urb[cpt], surf_ulb[cpt], surf_dlb[cpt]};
      Physical Volume(6000+cpt) = {vl[cpt-1]};
    EndIf

    If (cpt > 1 && cpt < N_DOM)
      Physical Surface(-(3+cpt)*1000) = {surf_drf[cpt-1], surf_urf[cpt-1],
        surf_ulf[cpt-1], surf_dlf[cpt-1], surf_drb[cpt-1], surf_urb[cpt-1],
        surf_ulb[cpt-1], surf_dlb[cpt-1]};
      Physical Surface((4+cpt)*1000) = {surf_drf[cpt], surf_urf[cpt], surf_ulf[cpt],
        surf_dlf[cpt], surf_drb[cpt], surf_urb[cpt], surf_ulb[cpt], surf_dlb[cpt]};
      Physical Volume(6000+cpt) = {vl[cpt-1]};
    EndIf

    Printf("Meshing sphere subdomain %g...", cpt-1);
    Save StrCat(MSH_NAME, Sprintf("%g.msh", cpt-1));

  EndFor
EndIf

// Physicals for the full domain
/*
Delete Physicals;
Physical Surface(-1) = {surf_drf[0], surf_urf[0], surf_ulf[0], surf_dlf[0],
                        surf_drb[0], surf_urb[0], surf_ulb[0], surf_dlb[0]};
Physical Surface(2) = {surf_drf[N_DOM], surf_urf[N_DOM], surf_ulf[N_DOM],
                       surf_dlf[N_DOM], surf_drb[N_DOM], surf_urb[N_DOM],
                       surf_ulb[N_DOM], surf_dlb[N_DOM]};
Physical Volume(6) = {vl[]};
Save "sphere_concentric.msh";
*/

BoundingBox {-R_EXT, R_EXT, -R_EXT, R_EXT, 0, 0};