/*********************************************************************
 *
 * This is a cylindrical mesh script for gmsh
 * the aim is to use it as mesh template for cylindrical reactor mesh generation
 *
 *********************************************************************/

u__1   = 0.001; //convert to meters
cl1    = u__1*10;   //mesh size at point ...

R2     = 125; // reactor radius
height = 200; // extrude in z direction
pz0    = 0;   // reactor elevation

r = 50;  // Note: 0 < r < R2 ("inner radius")


/* ** Structured mesh configuration  ** */
arcs_points          = 10 ;   //radial resolution
outer_radial_points  = 5 ;
no_z_cells           = 10;  //layers


/*********************************************************************
*********************************************************************/

// basic points (generators)
p_1 = newp;    Point(p_1) = {0,       0,      0,        cl1};  // center of reactor


// central geometry
p_2 = newp;    Point(p_2) = {u__1*r,  0,      u__1*pz0, cl1};  // arc points
p_3 = newp;    Point(p_3) = {0,       u__1*r, u__1*pz0, cl1};
p_4 = newp;    Point(p_4) = {-u__1*r, 0,      u__1*pz0, cl1};
p_5 = newp;    Point(p_5) = {0,      -u__1*r, u__1*pz0, cl1};

g = 1.5; // kinda gravity to influence arcs amplitude
p__2 = newp;    Point(p__2) = {u__1*r/g,   u__1*r/g,        u__1*pz0, cl1};  // outer points  influencing arcs amplitude
p__3 = newp;    Point(p__3) = {-u__1*r/g,   u__1*r/g,        u__1*pz0, cl1};  // arc points
p__4 = newp;    Point(p__4) = {-u__1*r/g,   -u__1*r/g,        u__1*pz0, cl1};  // arc points
p__5 = newp;    Point(p__5) = {u__1*r/g,   -u__1*r/g,        u__1*pz0, cl1};  // arc points

 // central arcs using BSpline 
lc_1=newl;    BSpline(lc_1) = {p_2,p__2,p_3};  
lc_2=newl;    BSpline(lc_2) = {p_3,p__3,p_4}; 
lc_3=newl;    BSpline(lc_3) = {p_4,p__4,p_5}; 
lc_4=newl;    BSpline(lc_4) = {p_5,p__5,p_2}; 

 
// central arcs using Circles
//lc_1=newl;    Circle(lc_1) = {p_2,p_1,p_3}; 
//lc_2=newl;    Circle(lc_2) = {p_3,p_1,p_4};
//lc_3=newl;    Circle(lc_3) = {p_4,p_1,p_5};
//lc_4=newl;    Circle(lc_4) = {p_5,p_1,p_2};


//Central surface
ll_1 = newll; 
Line Loop(ll_1)  = {p_1,p_2,p_3,p_4};

s1=news;
Ruled Surface(s1) = {ll_1};



/*****************************************************************************/
// Outer geometry
// circumference generators
p_22 = newp;    Point(p_22) = {u__1*R2,   0,        u__1*pz0, cl1};  // arc points
p_23 = newp;    Point(p_23) = {0,         u__1*R2,  u__1*pz0, cl1};
p_24 = newp;    Point(p_24) = {-u__1*R2,  0,        u__1*pz0, cl1};
p_25 = newp;    Point(p_25) = {0,         -u__1*R2, u__1*pz0, cl1};


//arcs (a circle is made of 4 arcs)
lc_21=newl;    Circle(lc_21) = {p_22,p_1,p_23}; 
lc_22=newl;    Circle(lc_22) = {p_23,p_1,p_24};
lc_23=newl;    Circle(lc_23) = {p_24,p_1,p_25};
lc_24=newl;    Circle(lc_24) = {p_25,p_1,p_22};

// lines to make the quadrants
lr_1=newl;    Line(lr_1) = {p_2, p_22};
lr_2=newl;    Line(lr_2) = {p_3, p_23};
lr_3=newl;    Line(lr_3) = {p_4, p_24};
lr_4=newl;    Line(lr_4) = {p_5, p_25};


// 1st quadrant
ll_q1 = newll;
Line Loop(ll_q1) = {lr_1, lc_21, -lr_2, -lc_1};
sq_1=news;
Plane Surface(sq_1) = {ll_q1};



// 2nd  quadrant
ll_q2 = newll;
//Line Loop(ll_q2) = {lr_1, lc_21, -lr_2, -lc_1};
Line Loop(ll_q2) = {lr_2, lc_22, -lr_3, -lc_2};
sq_2=news;
Plane Surface(sq_2) = {ll_q2};
//



// 3th quadrant
ll_q3 = newll;
Line Loop(ll_q3) = {lr_3, lc_23, -lr_4, -lc_3};
sq_3=news;
Plane Surface(sq_3) = {ll_q3};
//



// 4th quadrant
ll_q4 = newll;
Line Loop(ll_q4) = {lr_4, lc_24, -lr_1, -lc_4};
sq_4=news;
Plane Surface(sq_4) = {ll_q4};
//




//        Refining         //
//////////////////////////////////////////////////////////////
// outer radial refining
Transfinite Line {lr_1,lr_2,lr_3,lr_4}  = outer_radial_points;
Transfinite Line {lc_24,-lc_4} = arcs_points;
Transfinite Line {lc_23,-lc_3} = arcs_points;
Transfinite Line {lc_22,-lc_2} = arcs_points;
Transfinite Line {lc_21,-lc_1} = arcs_points;
Transfinite Line{lc_1,lc_2,lc_3,lc_4} =  arcs_points;
////////////////////////////////////////////////////////////



/****************************************************************
****************************************************************/
// joining al the surfaces


/* Make a list of a ring (annulus) of surfaces */
allParts[] = {s1,sq_1,sq_2,sq_3,sq_4};
 

/* Make surfaces to be meshed by transfinite algorithm */
Transfinite Surface {allParts[]};
 
/* The "Recombine Surface" command is issued in order to
 * create quadrilateral elements. */
Recombine Surface {allParts[]};


/****************************************************************
****************************************************************/
// Extruding circle to do a cylinder
/* Extrude geometry and quadrilateral mesh in Z direction */
//

//FELIX: I commented from here
/*
extrude_zdir[]=
Extrude{0,0,u__1*(height+pz0)} {
  Surface{allParts[]};
  Layers{no_z_cells};
  Recombine;
  };
  
  
num_volumes = 0;
num_surfs   = 0;
num_surfs_wall = 0;

    // 6 entries created by extrusion of a surface:
  //    [0] = created opposite surface
  //    [1] = created volume
  // [2..5] = created side surfaces
For i In {0:4}
    volumes[num_volumes]        = extrude_zdir[i*6+1]; num_volumes++;
	//Transfinite Volume{volumes[num_volumes-1]}={};
    surfs[num_surfs]            = extrude_zdir[i*6];   num_surfs++;
	//Surface {surfs[num_surfs-1]} In Volume {num_volumes-1};
	//Transfinite Surface {surfs[num_surfs-1]}={};
    walls_surfs[num_surfs_wall] = extrude_zdir[i*6+2]; num_surfs_wall++;
    walls_surfs[num_surfs_wall] = extrude_zdir[i*6+3]; num_surfs_wall++;
    walls_surfs[num_surfs_wall] = extrude_zdir[i*6+4]; num_surfs_wall++;
    walls_surfs[num_surfs_wall] = extrude_zdir[i*6+5]; num_surfs_wall++;
EndFor


withDilate = 1;

If(withDilate==1)
  //Dilate {{0, 0, u__1*(pz0)}, 1.2} { Surface{ surfs[] }; }
  /*FELIX: I changed this to get the surface at the same Z for other factors
  Not sure if that is the intended purpose/*
  Dilate {{0, 0, u__1*(height+pz0)}, 1.2} { Surface{ surfs[] }; } 
EndIf
*/ // FELIX: End of my commented section

// END of script

// FELIX: Copy and Dilate
Geometry.CopyMeshingMethod=1;
copy_dilate[]=
Dilate {{0, 0, u__1*(height+pz0)}, 1.2} { //Not sure about this. Check comment on L196
  Translate {0, 0, u__1*(height+pz0)} {
      Duplicata { Surface{ allParts[] }; }
  }
};

/* FELIX: This was an attempt to build a User-defined macro to create the volume connection
between the original and the dilated surface. I started having problem with duplicate entities.
 */
/*
surf_orig = 1;
surf_dila = 2;
Macro TransDilatedVolume
   b_lin_orig[]=Boundary{Surface{surf_orig};};
   b_lin_dila[]=Boundary{Surface{surf_dila};};
   For iL In {0:3}
      b_nod_orig[]=Boundary{Line{b_lin_orig[iL]};};
      b_nod_dila[]=Boundary{Line{b_lin_dila[iL]};};
      lin_1 = newl;
      Line(lin_1)={b_nod_orig[0],b_nod_dila[0]};
      //lin_2 = newl;
      //Line(lin_2)={b_nod_orig[1],b_nod_dila[1]};
      Transfinite Line{lin_1}=no_z_cells+1 Using Progression 1;
      //Printf("Line Loop: %g %g %g %g",b_lin_orig[iL],lin_1,b_lin_dila[iL],-lin_2);
      //newsurf = newll;
      //Line Loop(newsurf) = {b_lin_orig[iL],lin_1,b_lin_dila[iL],-lin_2};
      //Ruled Surface(newsurf)={newsurf};
   EndFor
Return

//Call of the Macro
surf_orig=sq_1;
surf_dila=copy_dilate[0];
Call TransDilatedVolume;

//Printing info
Printf("WHICH Surface sq1? : %g",sq_1);
Printf("WHICH Lines? :");
tmp[]=Boundary{Surface{sq_1};};
For i In {0:3}
   Printf("%g",tmp[i]);
EndFor

For i In {0:4}
   Printf("Surface: %g",copy_dilate[i]);
   b_lin[]=Boundary{Surface{copy_dilate[i]};};
   For j In {0:3}
      Printf("   Line %g",b_lin[j]);
      b_nod[]=Boundary{Line{b_lin[j]};};
      For k In {0:1}
         Printf("      Vertex %g",b_nod[k]);
      EndFor
   EndFor
EndFor
*/

//FELIX: This section was added by the GUI when I did the stuff I said in my mail
/*
Line(46) = {11, 40};
Line(47) = {12, 58};
Line(48) = {4, 21};
Line(49) = {3, 16};

Line Loop(50) = {47, 36, -48, 13};
Line Loop(51) = {46, 31, -49, 12};
Plane Surface(52) = {50};
Plane Surface(53) = {51};
Line Loop(54) = {25, -48, -2, 49};
Ruled Surface(55) = {54};
Line Loop(56) = {35, -47, -8, 46};
Ruled Surface(57) = {56};

Transfinite Line {47, 46, 48, 49} = 10 Using Progression 1;
Transfinite Surface {52};
Transfinite Surface {53};
Transfinite Surface {57};
Transfinite Surface {55};
Recombine Surface {52, 53, 57, 55};
Surface Loop(58) = {57, 33, 52, 55, 53, 18};
Volume(59) = {58};
Transfinite Volume {59};
*/
Mesh.RecombineAll = 1;