// to preserve transfinite in geometry duplications
Geometry.CopyMeshingMethod = 1 ;

Function Quadrangle

  /*
    * X, Y, Z, lc
    * point   dX    dY
    * 1       0     0
    * 2       dX2   dY2
    * 3       dX3   dY3
    * 4       dX4   dY4
    *
    * Transfinite:
    * Horizontal: N_H
    * Vertical: N_V
    */
  p1 = newp; Point(p1) = {X,Y,Z,lc};
  p2 = newp; Point(p2) = {X+dX2,Y+dY2,Z,lc};
  p3 = newp; Point(p3) = {X+dX3,Y+dY3,Z,lc};
  p4 = newp; Point(p4) = {X+dX4,Y+dY4,Z,lc};
  
  l1= newl; Line(l1) = {p1, p2} ;
  l2= newl; Line(l2) = {p2, p3} ;
  l3= newl; Line(l3) = {p3, p4} ;
  l4= newl; Line(l4) = {p4, p1} ;
  
  Transfinite Line {l1, l3} = N_H;
  Transfinite Line {l2, l4} = N_V;
  
  ll1 = newll; Line Loop(ll1) = {l1, l2, l3, l4};
  
  pl1 = news;
  
  Plane Surface (pl1) = {ll1};
  
  Transfinite Surface {pl1};
  
  Recombine Surface {pl1};

Return


Function Rectangle

  dX3=dX2;
  dX4=0;
  dY2=0;
  dY4=dY3;
  Call Quadrangle ;

Return


Function Square

  dY3=dX2;
  N_V = N_H;
  Call Rectangle ;
  
Return


Function drawCircle

  center = newp; Point(center) = {X,Y,Z,lc};
  helpVar = D/2 * 0.70710678/*sin (45)*/; 
  p1 = newp; Point(p1) = {X + helpVar, Y ,Z + helpVar,lc};
  p2 = newp; Point(p2) = {X - helpVar, Y ,Z + helpVar,lc};
  p3 = newp; Point(p3) = {X - helpVar, Y ,Z - helpVar,lc};
  p4 = newp; Point(p4) = {X + helpVar, Y ,Z - helpVar,lc};
  
  c1= newc; Circle(c1) = {p1, center, p2} ;
  c2= newc; Circle(c2) = {p2, center, p3} ;
  c3= newc; Circle(c3) = {p3, center, p4} ;
  c4= newc; Circle(c4) = {p4, center, p1} ;

  Transfinite Line {c1, c2, c3, c4} = N_H Using Bump 0.5;

  ll1 = newll; Line Loop(ll1) = {c1, c2, c3, c4};
  
  pl1 = news; Plane Surface (pl1) = {ll1};

  Transfinite Surface {pl1};
  
  Recombine Surface {pl1};

Return


Function drawHalfCircle

  center = newp; Point(center) = {X,Y,Z,lc};
  helpVar = D/2 * 0.70710678/*sin (45)*/; 
  p1 = newp; Point(p1) = {X + helpVar, Y ,Z + helpVar,lc};
  p2 = newp; Point(p2) = {X - helpVar, Y ,Z + helpVar,lc};
  p3 = newp; Point(p3) = {X - D/2, Y ,Z,lc};
  p4 = newp; Point(p4) = {X + D/2, Y ,Z,lc};
  
  c1= newc; Circle(c1) = {p1, center, p2} ;
  c2= newc; Circle(c2) = {p2, center, p3} ;
  l3= newl; Line(l3) = {p3, p4} ;
  c4= newc; Circle(c4) = {p4, center, p1} ;

  Transfinite Line {c1, c2, l3, c4} = N_H Using Bump 0.5;

  ll1 = newll; Line Loop(ll1) = {c1, c2, l3, c4};
  
  pl1 = news; Plane Surface (pl1) = {ll1};

  Transfinite Surface {pl1};
  
  Recombine Surface {pl1};

Return





Function ExtrudeFunc 
 
  myExt[]=Extrude {ExX, ExY, ExZ}{
      Surface{pl1};
      Layers{ExLayers};
      Recombine;
  };

Return

Function ExtrudeFuncNoTrans 
 
  myExt[]=Extrude {ExX, ExY, ExZ}{
      Surface{pl1};
      Recombine;
  };
  Recombine;

Return

Function ExtrudeRotFunc 

  myExt[]=Extrude {{0, 0, 1}, {X,Y,Z}, Pi}{
      Surface{pl1};
      Layers{ExLayers};
      Recombine;
  };

Return


Function Cuboid
  Call Quadrangle;
  Call ExtrudeFunc;
Return

Function RectCuboid
  Call Rectangle;
  Call ExtrudeFunc;
Return

Function Cube
  Call Square;
  ExLayers = N_H;
  ExZ=dX2;
  Call ExtrudeFunc;
Return

Function Cylinder
  Call drawCircle;
  Call ExtrudeFunc;
Return

Function halfCylinder
  Call drawHalfCircle;
  Call ExtrudeFunc;
Return


Function pileHouseSketch
  center = newp; Point(center) = {X,Y,Z,lc};
  helpVar = D/2 * 0.70710678/*sin (45)*/; 
  p1 = newp; Point(p1) = {X + helpVar, Y ,Z + helpVar,lc};
  p2 = newp; Point(p2) = {X - helpVar, Y ,Z + helpVar,lc};
  c1= newc; Circle(c1) = {p1, center, p2} ;
  outerLength = outerLengthDiamRatio * D;
helpVar2 = outerLength/2 ;//* 0.70710678/*sin (45)*/; 
  p3 = newp; Point(p3) = {X + helpVar2, Y ,Z + helpVar2,lc};
  p4 = newp; Point(p4) = {X - helpVar2, Y ,Z + helpVar2,lc};
  l1= newl; Line (l1) = {p2, p4};
  l2= newl; Line (l2) = {p4, p3};
  l3 = newl; Line (l3) = {p3, p1};
  Transfinite Line {c1} = N_H Using Bump 0.5;
  Transfinite Line {l2} = N_H;
  Transfinite Line {l1,l3} = N_V;
  ll1 = newll; Line Loop(ll1) = {c1, l1, l2, l3};
    
  pl1 = news; Plane Surface (pl1) = {ll1};

  Transfinite Surface {pl1};
  
  Recombine Surface {pl1};
Return

Function halfPileHouseSketch
  center = newp; Point(center) = {X,Y,Z,lc};
  helpVar = D/2 * 0.70710678/*sin (45)*/; 
  p1 = newp; Point(p1) = {X + helpVar, Y ,Z + helpVar,lc};
  p2 = newp; Point(p2) = {X + D/2, Y ,Z,lc};
  c1= newc; Circle(c1) = {p1, center, p2} ;
  outerLength = outerLengthDiamRatio * D;
helpVar2 = outerLength/2 ;//* 0.70710678/*sin (45)*/; 
  p3 = newp; Point(p3) = {X + helpVar2, Y ,Z + helpVar2,lc};
p4 = newp; Point(p4) = {X + helpVar2/*0.70710678*/, Y ,Z,lc};
  l1= newl; Line (l1) = {p2, p4};
  l2= newl; Line (l2) = {p4, p3};
  l3 = newl; Line (l3) = {p3, p1};
  Transfinite Line {c1} = N_H Using Bump 0.5;
  Transfinite Line {l2} = N_H;
  Transfinite Line {l1,l3} = N_V;
   ll1 = newll; Line Loop(ll1) = {c1, l1, l2, l3};
  
  pl1 = news; Plane Surface (pl1) = {ll1};

  Transfinite Surface {pl1};
  
  Recombine Surface {pl1};
Return

Function otherHalfPileHouseSketch
  center = newp; Point(center) = {X,Y,Z,lc};
  helpVar = D/2 * 0.70710678/*sin (45)*/; 
  p1 = newp; Point(p1) = {X - helpVar, Y ,Z + helpVar,lc};
  p2 = newp; Point(p2) = {X - D/2, Y ,Z,lc};
  c1= newc; Circle(c1) = {p1, center, p2} ;
  outerLength = outerLengthDiamRatio * D;
helpVar2 = outerLength/2 ;//* 0.70710678/*sin (45)*/; 
  p3 = newp; Point(p3) = {X - helpVar2, Y ,Z + helpVar2,lc};
p4 = newp; Point(p4) = {X - helpVar2/*0.70710678*/, Y ,Z,lc};
  l1= newl; Line (l1) = {p2, p4};
  l2= newl; Line (l2) = {p4, p3};
  l3 = newl; Line (l3) = {p3, p1};
  Transfinite Line {c1} = N_H Using Bump 0.5;
  Transfinite Line {l2} = N_H;
  Transfinite Line {l1,l3} = N_V;
   ll1 = newll; Line Loop(ll1) = {c1, l1, l2, l3};
  
  pl1 = news; Plane Surface (pl1) = {ll1};

  Transfinite Surface {pl1};
  
  Recombine Surface {pl1};
Return


/*Function Sphere
    Call drawCircle;
    Call ExtrudeRotFunc;
    Return*/

Function pileHouse
  Call pileHouseSketch;
  Call ExtrudeFunc;
  myPi=3.14159265359;
  pileHouse1 = myExt[1]; 
//Symmetry { 0.0,0.0,1.0,0.0 }{Duplicata{Surface{pl1};}}
  rot1[] = Rotate{{0,1,0},{X,Y,Z},myPi/2}{Duplicata{Surface{pl1};}};
  rot2[] = Rotate{{0,1,0},{X,Y,Z},myPi}{Duplicata{Surface{pl1};}};
  rot3[] = Rotate{{0,1,0},{X,Y,Z},-myPi/2}{Duplicata{Surface{pl1};}};
  pl1 = rot1[0];
  Call ExtrudeFunc;
  pileHouse2 = myExt[1];
  pl1 = rot2[0];
  Call ExtrudeFunc;
  pileHouse3 = myExt[1];
  pl1 = rot3[0];
  Call ExtrudeFunc;
  pileHouse4 = myExt[1];
Return


Function halfPileHouse
  Call pileHouseSketch;
  Call ExtrudeFunc;
  myPi=3.14159265359;
  pileHouse1 = myExt[1]; 
  Call halfPileHouseSketch;
  Call ExtrudeFunc;
  pileHouse2 = myExt[1];
  Call otherHalfPileHouseSketch;
  Call ExtrudeFunc;
  pileHouse3 = myExt[1];
Return



Function HorzQuadrangle

  /*
    * X, Y, Z, lc
    * point   dX    dZ
    * 1       0     0
    * 2       dX2   dZ2
    * 3       dX3   dZ3
    * 4       dX4   dZ4
    *
    * Transfinite:
    * Horizontal: N_H (X)
    * Vertical: N_V (Z)
    */
  p1 = newp; Point(p1) = {X,Y,Z,lc};
  p2 = newp; Point(p2) = {X+dX2,Y,Z+dZ2,lc};
  p3 = newp; Point(p3) = {X+dX3,Y,Z+dZ3,lc};
  p4 = newp; Point(p4) = {X+dX4,Y,Z+dZ4,lc};
  
  l1= newl; Line(l1) = {p1, p2} ;
  l2= newl; Line(l2) = {p2, p3} ;
  l3= newl; Line(l3) = {p3, p4} ;
  l4= newl; Line(l4) = {p4, p1} ;
  
  Transfinite Line {l1, l3} = N_H;
  Transfinite Line {l2, l4} = N_V;
  
  ll1 = newll; Line Loop(ll1) = {l1, l2, l3, l4};
  
  pl1 = news;
  
  Plane Surface (pl1) = {ll1};
  
  Transfinite Surface {pl1};
  
  Recombine Surface {pl1};

Return



Function HorzQuadrangleNoTrans

  /*
    * X, Y, Z, lc
    * point   dX    dZ
    * 1       0     0
    * 2       dX2   dZ2
    * 3       dX3   dZ3
    * 4       dX4   dZ4
    *
    * Transfinite:
    * Horizontal: N_H (X)
    * Vertical: N_V (Z)
    */
  p1 = newp; Point(p1) = {X,Y,Z,lc};
  p2 = newp; Point(p2) = {X+dX2,Y,Z+dZ2,lc};
  p3 = newp; Point(p3) = {X+dX3,Y,Z+dZ3,lc};
  p4 = newp; Point(p4) = {X+dX4,Y,Z+dZ4,lc};
  
  l1= newl; Line(l1) = {p1, p2} ;
  l2= newl; Line(l2) = {p2, p3} ;
  l3= newl; Line(l3) = {p3, p4} ;
  l4= newl; Line(l4) = {p4, p1} ;
  
//  Transfinite Line {l1, l3} = N_H;
//  Transfinite Line {l2, l4} = N_V;
  
  ll1 = newll; Line Loop(ll1) = {l1, l2, l3, l4};
  
  pl1 = news;
  
  Plane Surface (pl1) = {ll1};
  
//  Transfinite Surface {pl1};
  
  Recombine Surface {pl1};

Return


Function HorzRectangle

  dX3=dX2;
  dX4=0;
  dZ2=0;
  dZ4=dZ3;
  Call HorzQuadrangle ;

Return

Function HorzRectangleNoTrans

  dX3=dX2;
  dX4=0;
  dZ2=0;
  dZ4=dZ3;
  Call HorzQuadrangleNoTrans ;

Return


Function HorzSquare

  dZ3=dX2;
  N_V = N_H;
  Call HorzRectangle ;
  
Return


Function HorzCuboid
  Call HorzQuadrangle;
  Call ExtrudeFunc;
Return

Function HorzRectCuboid
  Call HorzRectangle;
  Call ExtrudeFunc;
Return

Function HorzRectCuboidNoTrans
  Call HorzRectangleNoTrans;
//Call ExtrudeFunc;
  Call ExtrudeFuncNoTrans;
Return


Function HorzCube
  Call HorzSquare;
  ExLayers = N_H;
  ExY=dX2;
  Call ExtrudeFunc;
Return


/////
lc = 0.1;
/////
X = 0; 
/////
Y = 0;
/////
Z = 0;
/////
dX2 = 1;
/////
dX3 = 0.9;
/////
dX4= 0.4;
/////
dY2= 0.05;
/////
dY3= 0.2;
/////
dY4= 0.18;

/////
N_V=25;
/////
N_H=20;

//Call quadrangle ; 

//Call Rectangle ;

//Call Square ;

/////
D=0.05;
/////
N_H=31;
//Call drawCircle ;

/////
ExX = 0;
/////
ExY = 1;
/////
ExZ = 0;

/////
ExLayers = 30;
/////
outerLengthDiamRatio = 2;// 2D

//Call Cube;

//Call Cuboid;

//Call RectCuboid;

//Call  

//Call pileHouse; 
//Call halfPileHouse;
//Call Cylinder;
//Call halfCylinder;
//Call Sphere;

// Automatically Assign Physical Volume to the Extrude
// Physical Volume ("myVolName") = {myExt[1]}; 
// myExt[0] is the ID of the opposite new surface

// Physical Surfaces are to be defined manually

Mesh.Smoothing = 100;
