//Geometry
forcoolant = 0;
forfuel = 0;
forGraphit = 1;

R = 0.2;
r = 0.1;
pitch = .4;
gap = 0;

height = 0.1;
coreNodes = 1;
rows = 5;

edge = 2;
halfEdge = edge/2;
flatToFlat = edge * Sqrt(3);
halfFlatToFlat = flatToFlat/2;

lc = 0.1;
// Coolant Position
Cposxn = {};
Cposyn = {};

Cposxn += {2};
Cposyn += {0};

Cposxn += {3};
Cposyn += {-1};

// Cposxn += {4};
// Cposyn += {-1};

// Cposxn += {4};
// Cposyn += {1};

Cposx = {};
Cposy = {};
CposX = {};
CposY = {};
For i In {0 : (#Cposxn[]-1)}
  Cposx += {Cposxn[i]*pitch};
  Cposy += {Cposyn[i]*pitch};
  CposX += {Cposx[i] * Sin(Pi/3)};
  CposY += {Cposx[i] * Cos(Pi/3) + Cposy[i]};
  Printf("CposX %g and CposY %g", CposX[i], CposY[i]);

      For kk In {1:5}
      CposX0 = 0.5*CposX[#CposX[]-1] - Sin(Pi/3)*CposY[#CposX[]-1];
      CposY0 = Sin(Pi/3)*CposX[#CposX[]-1] + 0.5*CposY[#CposX[]-1];
      CposX += CposX0;
      CposY += CposY0;
      Printf("CposX %g and CposY %g", CposX[i*5+kk], CposY[i*5+kk]);

   	  EndFor  
EndFor



// Circle
Point(1) = {0,0,-1,0.5};

Point(2) = {r,0,-1,0.5};
Point(3) = {0,r,-1,0.5};
Point(4) = {-r,0,-1,0.5};
Point(5) = {0,-r,-1,0.5};

Circle(1) = {2,1,3};
Circle(2) = {3,1,4};
Circle(3) = {4,1,5};
Circle(4) = {5,1,2};

Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};


Point(6) = {R,0,-1,0.5};
Point(7) = {0,R,-1,0.5};
Point(8) = {-R,0,-1,0.5};
Point(9) = {0,-R,-1,0.5};

Circle(7) = {6,1,7};
Circle(8) = {7,1,8};
Circle(9) = {8,1,9};
Circle(10) = {9,1,6};

Line Loop(11) = {7,8,9,10};
Plane Surface(12) = {11};

Point(11) = {0, -edge,  -1, lc} ;
Point(12) = {halfFlatToFlat, -halfEdge, -1, lc} ;
Point(13) = {halfFlatToFlat, halfEdge,  -1, lc} ;
Point(14) = {0, edge,  -1, lc} ;
Point(15) = {-halfFlatToFlat, halfEdge, -1, lc} ;
Point(16) = {-halfFlatToFlat, -halfEdge,  -1, lc} ;

Line(17) = {11,12} ;
Line(18) = {12,13} ;
Line(19) = {13,14} ;
Line(20) = {14,15} ;
Line(21) = {15,16} ;
Line(22) = {16,11} ;

Line Loop(23) = {17,18,19,20,21,22} ;

Plane Surface(24) = {23} ;
surface[] = Translate {0, 0, 1} { Duplicata{ Surface{24}; } };
graphit[] = Extrude {0, 0, height} { 
Surface{surface[]};  Layers{coreNodes} ; Recombine;
};
GraphitVolumes[] = graphit[1];
GraphitTop[] = graphit[0];
GraphitBtm[] = surface[];
GraphitWall[] = graphit[2];
GraphitWall[] += graphit[3];
GraphitWall[] += graphit[4];
GraphitWall[] += graphit[5];
GraphitWall[] += graphit[6];
GraphitWall[] += graphit[7];


// Declare
CoolantVolumes[] = {};
pinVolumes[] = {};
CoolantTop[] = {};
CoolantBtm[] = {};
CoolantWall[] = {};
pinTop[] = {};
pinBtm[] = {};
pinWall[] = {};


// Function

Function insertPin
  // Printf("x %g and y %g", x, y);
  // Printf("X %g and Y %g", X, Y);
  surface[] = Translate {X, Y, 1} { Duplicata{ Surface{6}; } };
  // Printf("surface 1%g", surface[]);
  pin[] = Extrude {0, 0, height} { 
  Surface{surface[0]};  Layers{coreNodes} ; Recombine;
  };

  
	// If(gap != 0.0)
	// 	surfaceGap[] = Translate {X, Y, 1+lowReflHeight+lowFGPHeight+diagridHeight+coreHeight+upperFGPHeight} { Duplicata{ Surface{2}; } };
	// 	assemblyGap[] = Extrude {0, 0, upperReflHeight} { 
	// 	Surface{surfaceGap[0]};  Layers{upperReflNodes} ; Recombine;
	// 	};
	// 	assemblyGapVolumes += assemblyGap[1];
	// EndIf

 
  isCoolant = 0 ; 
  For l In {0 : (#CposX[] - 1)}
    // Printf("diffX2 %g and diffY2 %g", (X-CposX[l])*(X-CposX[l]), (Y-CposY[l])*(Y-CposY[l]));
    If(((X-CposX[l])*(X-CposX[l]))<0.01 && ((Y-CposY[l])*(Y-CposY[l]))<0.01)
      isCoolant = 1;
      // Printf("coolantX %g and coolantY %g", X, Y);
    EndIf
  EndFor
  
  If(isCoolant)
    CoolantVolumes += pin[1];
    CoolantTop += pin[0];
    Printf("CoolantTop +%g", pin[0]);
    CoolantBtm += surface[0];
    Printf("CoolantBtm +%g", surface[]);
    CoolantWall += pin[2];
    CoolantWall += pin[3];
    CoolantWall += pin[4];
    CoolantWall += pin[5];

  EndIf
  If(!isCoolant)
    pinVolumes += pin[1];
    pinTop += pin[0];
    pinBtm += surface[0];
    pinWall+= pin[2];
	pinWall+= pin[3];
	pinWall+= pin[4];
	pinWall+= pin[5];
  EndIf

Return

For j In {1 : (rows-1)}
  x = 0;
  y = j*pitch;
  X = x * Sin(Pi/3);
  Y = x * Cos(Pi/3) + y;
  Call insertPin; 
  If (j>0)
  For k In {1:5}
    X0 = 0.5*X - Sin(Pi/3)*Y;
    Y0 = Sin(Pi/3)*X + 0.5*Y;
    X = X0;
    Y = Y0;
    Call insertPin;      
  EndFor
  
  For i In {1:(j-1)}
    x = i*pitch;
    y -= pitch;
    X = x * Sin(Pi/3);
    Y = x * Cos(Pi/3) + y;
    Call insertPin;
    For k In {1:5}
      X0 = 0.5*X - Sin(Pi/3)*Y;
      Y0 = Sin(Pi/3)*X + 0.5*Y;
      X = X0;
      Y = Y0;
      Call insertPin;      
    EndFor
  EndFor
  EndIf
EndFor

X=0;
Y=0;
surface[] = Translate {X, Y, 1} { Duplicata{ Surface{12}; } };
pin[] = Extrude {0, 0, height} { 
Surface{surface[0]};  Layers{coreNodes} ; Recombine;
};
CoolantVolumes += pin[1];
CoolantTop += pin[0];
CoolantBtm += surface[0];
CoolantWall += pin[2];
CoolantWall += pin[3];
CoolantWall += pin[4];
CoolantWall += pin[5];

Coherence;

If(forcoolant)
Physical Volume("coolant") = {CoolantVolumes[]} ;
Physical Surface("frontAndBack") = {CoolantBtm[], CoolantTop[]} ;
// Physical Surface("bottom") = {CoolantBtm[]} ;
Physical Surface("walls") = {CoolantWall[]} ;
// Physical Surface("top") = {CoolantTop[]} ;
EndIf
If(forfuel)
Physical Volume("fuel") = {pinVolumes[]} ;
Physical Surface("frontAndBack") = {pinBtm[], pinTop[]} ;
// Physical Surface("bottom") = {pinBtm[]} ;
Physical Surface("walls") = {pinWall[]} ;
// Physical Surface("top") = {pinTop[]} ;
EndIf
If(forGraphit)
Physical Volume("Graphit") = {GraphitVolumes[], -pinVolumes[], -CoolantVolumes[]} ;
Physical Surface("frontAndBack") = {GraphitBtm[], -pinBtm[], -CoolantBtm[], GraphitTop[], -pinTop[], -CoolantTop[]} ;
// Physical Surface("bottom") = {GraphitBtm, -pinBtm, -CoolantBtm} ;
Physical Surface("walls") = {GraphitWall[], -pinWall[], -CoolantWall[]} ;
// Physical Surface("top") = {GraphitTop, -pinTop, -CoolantTop};
EndIf
 
Recombine;
Delete {Surface{6}; Line{1};  Line{2}; Line{3}; Line{4}; Point{1};  Point{2}; Point{3}; Point{4}; Point{5};};
Delete {Surface{12}; Line{8};  Line{9}; Line{10}; Line{7}; Point{1}; Point{6};  Point{7}; Point{8}; Point{9};};
Coherence;
Coherence;
