Merge "circleLine1.mesh";
lineBnd[0] = 10;

Merge "circleLine2.mesh";
lineBnd[1] = 20;

Merge "object.stl";
objSurf[] = Surface "*";

Coherence Mesh;


////////////////////////////////////////////////////////////////////////////////

Lcar1 = 0.1;
Lcar3 = 2.0;
Lcar5 = 5.0;

Mesh.Algorithm = 1;
Mesh.CharacteristicLengthExtendFromBoundary = 1;
Mesh.CharacteristicLengthMin = Lcar1*2/3;
Mesh.CharacteristicLengthMax = Lcar5;
Mesh.Optimize = 1;
Mesh.OptimizeNetgen = 1;

////////////////////////////////////////////////////////////////////////////////
// Create inlet surface

ll1 = newll;
Curve Loop(ll1) = {lineBnd[0]};
sInlet = news;
Plane Surface(sInlet) = {ll1};
Printf("Inlet: %g", sInlet);

////////////////////////////////////////////////////////////////////////////////
// Create a square with a circular hole inside

xx[] = {+50.0, -50.0, -50.0, +50.0};
yy[] = {+50.0, +50.0, -50.0, -50.0};

pp2[0] = newp;
Point(pp2[0]) = {xx[0], yy[0], 0, Lcar5};
pp2[1] = newp;
Point(pp2[1]) = {xx[1], yy[1], 0, Lcar5};
pp2[2] = newp;
Point(pp2[2]) = {xx[2], yy[2], 0, Lcar5};
pp2[3] = newp;
Point(pp2[3]) = {xx[3], yy[3], 0, Lcar5};

For ii In {0:2}
  rectT[ii] = newl;
  Line(rectT[ii]) = {pp2[ii],pp2[ii+1]};
EndFor
rectT[3] = newl;
Line(rectT[3]) = {pp2[3],pp2[0]};

rectLoopTOPwithHole=newll;
Curve Loop(rectLoopTOPwithHole)={rectT[],-lineBnd[1]};
TopBoxHole=news;
Plane Surface(TopBoxHole) = {rectLoopTOPwithHole};
Printf("TopExt: %g", TopBoxHole);

////////////////////////////////////////////////////////////////////////////////
// Create bottom and lateral surfaces of the box

zz[] = {-30.0 , -30.0 , -30.0 , -30.0};

pp2[4] = newp;
Point(pp2[4]) = {xx[0], yy[0], zz[0], Lcar5};
pp2[5] = newp;
Point(pp2[5]) = {xx[1], yy[1], zz[1], Lcar5};
pp2[6] = newp;
Point(pp2[6]) = {xx[2], yy[2], zz[2], Lcar5};
pp2[7] = newp;
Point(pp2[7]) = {xx[3], yy[3], zz[3], Lcar5};

For ii In {0:2}
  rectB[ii] = newl;
  Line(rectB[ii]) = {pp2[ii+4],pp2[ii+5]};
EndFor
rectB[3] = newl;
Line(rectB[3]) = {pp2[7],pp2[4]};

For ii In {0:3}
  rectL[ii] = newl;
  Line(rectL[ii]) = {pp2[ii],pp2[ii+4]};
EndFor

rectLoopBOTTOM=newll;
Curve Loop(rectLoopBOTTOM)={rectB[]};
BottomBox=news;
Plane Surface(BottomBox) = {rectLoopBOTTOM};
Printf("Bottom: %g", BottomBox);

For ii In{0:2}
  rectLoopLateral[ii]=newll;
  Curve Loop(rectLoopLateral[ii])={rectT[ii],rectL[ii+1],-rectB[ii],-rectL[ii]};
  LateralBox[ii]=news;
  Plane Surface(LateralBox[ii]) = {rectLoopLateral[ii]};
  Printf("Lateral[%g]: %g", ii, LateralBox[ii]);
EndFor
rectLoopLateral[3]=newll;
Curve Loop(rectLoopLateral[3])={rectT[3],rectL[0],-rectB[3],-rectL[3]};
LateralBox[3]=news;
Plane Surface(LateralBox[3]) = {rectLoopLateral[3]};
Printf("Lateral[3]: %g", LateralBox[3]);

////////////////////////////////////////////////////////////////////////////////
// Create Volume

sl=newsl;
Surface Loop(sl) = {sInlet, objSurf[], TopBoxHole, LateralBox[], BottomBox};
vol=newv;
Volume(vol) = {sl};
//ReorientMesh Volume {vol};
