SetFactory("OpenCASCADE");
Geometry.Surfaces = 0;
Geometry.MatchMeshTolerance = 1.0e-8;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthFromCurvature = 0;
Mesh.CharacteristicLengthExtendFromBoundary = 0;

DefineConstant[alpha = {Pi/8.0, Name "alpha"}];
DefineConstant[w_over_D = {6.0, Name "w_over_D"}];

D = 30.0e-3;
R = D/2.0;
Linlet = 5.0*D;
Loutlet = 15.0*D;
L = Loutlet;
w = w_over_D*D;
dist = w/2.0;

// ------------------------------
// Create the cylinders
// ------------------------------

disk1 = news; Disk(disk1) = {0.0,0.0,0.0,R};
Rotate {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, (Pi/2.0)} { Surface{disk1}; }
Rotate {{0.0, 1.0, 0.0}, {0.0, 0.0, 0.0}, (Pi)} { Duplicata{Surface{disk1};} }
disk = news;
BooleanUnion(disk) = { Surface{disk1}; Delete; }{ Surface{disk1+1}; Delete; };

Translate {0.0,-w,0.0}{Duplicata {Surface{disk};}}
Translate {0.0,w,0.0}{Duplicata {Surface{disk};}}

ll0 = newll;
e() = Unique(Abs(Boundary{ Surface{disk}; }));
Curve Loop(ll0) = e();

ll1 = newll;
e() = Unique(Abs(Boundary{ Surface{disk+1}; }));
Curve Loop(ll1) = e();

ll2 = newll;
e() = Unique(Abs(Boundary{ Surface{disk+2}; }));
Curve Loop(ll2) = e();

cyl1 = newv;
ThruSections(cyl1) = {ll1,ll0,ll2};

// Rotate them at the alpha angle
Rotate {{0.0, 0.0, 1.0}, {0.0, 0.0, 0.0}, -alpha} {Duplicata{Volume{cyl1};} }
Rotate {{0.0, 0.0, 1.0}, {0.0, 0.0, 0.0}, alpha} {Volume{cyl1}; }

cyl = newv; BooleanUnion(cyl) = { Volume{cyl1+1}; Delete; }{ Volume{cyl1}; Delete; };

// ------------------------------
// Create diamond-shape fluid domain
// ------------------------------

X = dist * Sin(alpha);
Y = dist * Cos(alpha);

p0 = newp; Point(p0) = {-2.0*X,0.0,0.0};
p1 = newp; Point(p1) = {0.0,-2.0*Y,0.0};
p2 = newp; Point(p2) = {2.0*X,0.0,0.0};
p3 = newp; Point(p3) = {0.0,2.0*Y,0.0};

l0 =newl; Line(l0) = {p0,p1};
l1 =newl; Line(l1) = {p1,p2};
l2 =newl; Line(l2) = {p2,p3};
l3 =newl; Line(l3) = {p3,p0};

cl4 = newll; Curve Loop(cl4) = {l0:l3};

surf = news; Plane Surface(surf) = {cl4};

Translate{0.0,0.0,-Linlet}{Surface{surf};}
Extrude {0.0,0.0,Linlet+Loutlet}{Surface{surf};}

domain = newv;BooleanDifference(domain) = { Volume{cyl+1};Delete;}{ Volume{cyl}; Delete;};

Rotate {{0.0, 0.0, 1.0}, {0.0, 0.0, 0.0}, -alpha} {Volume{domain}; }

// ------------------------------
// Mesh size fields
// ------------------------------

//lc = Pi*D/30.0;
lc = Pi*D/15.0;

Field[1] = Box;
Field[1].VIn = 5.0*lc;
Field[1].VOut = 1000.0*lc;
Field[1].ZMin = -L;
Field[1].ZMax = L;
Field[1].XMin = -L;
Field[1].XMax = L;
Field[1].YMin = -L;
Field[1].YMax = L;

rect = news; Rectangle(rect)={-2.0*L,-2.0*L,0.0,4.0*L,4.0*L};
Field[2] = Attractor;
Field[2].NNodesByEdge = 100;
Field[2].FacesList = {rect};

//// We then define a Threshold field, which uses the return value of the
//// Attractor Field[2] in order to define a simple change in element size around
//// the attractors 
////
//// LcMax -                         /------------------
////                               /
////                             /
////                           /
//// LcMin -o----------------/
////        |                |       |
////     Attractor       DistMin   DistMax

Field[3] = Threshold;
Field[3].IField = 2;
Field[3].LcMin = 0.75*lc;
Field[3].LcMax = 3.0*lc;
Field[3].DistMin = 5.0*R;
Field[3].DistMax = 15.0*R;
//Field[3].Sigmoid = 1;
Field[3].StopAtDistMax = 1;

Field[4] = Min;
Field[4].FieldsList = {1,3};

Field[5] = Restrict;
Field[5].RegionsList = {domain};

s() = Unique(Abs(Boundary{Volume{domain};}));
Field[5].FacesList = {s()};

e() = Unique(Abs(Boundary{Surface{s()};}));
Field[5].EdgesList = {e()};

Field[5].IField = 4;

Background Field = 5;

Mesh.Algorithm = 6;
Mesh.Algorithm3D = 1;

Mesh.Optimize = 1;
Mesh.OptimizeNetgen = 1;

Mesh.CharacteristicLengthMax = 5.0*lc ;
Mesh.CharacteristicLengthMin = lc/1.5 ;

// ------------------------------
// Boundary conditions
// ------------------------------

Physical Volume("volume") = {domain};

Physical Surface("inlet") = {8};
Physical Surface("outlet") = {9};

Physical Surface("Periox1") = {6};
Physical Surface("Periox2") = {13};
Physical Surface("Perioy1") = {10};
Physical Surface("Perioy2") = {7};

Physical Surface("w_c_red") =    {17,16,14,15};
Physical Surface("w_c_blue") =   {11,12,18,19};

// ------------------------------
// Apply periodicities on specific boundaries
// ------------------------------
X = w * Cos(Pi/2.0 - 2.0*alpha);
Y = w * Sin(Pi/2.0 - 2.0*alpha);

Periodic Surface{7} = {10} Translate{0.0,w,0.0};
Periodic Surface{13} = {6} Translate{X,Y,0.0};
