Radius = 0.5;

lc = 4;
ht = 0.5;

Point(1) = { 0.0 , 0.0 , 0.0 ,    ht}; 		
Point(2) = { Radius , 0.0 , 0.0 , ht};
Point(3) = { 0.0 , Radius , 0.0 , ht};
Point(4) = { -Radius ,0.0 , 0.0 , ht};
Point(5) = {0.0 , -Radius , 0.0 , ht};

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

Point(6) = { -20 , -20 , 0.0 , lc}; 
Point(7) = { -20 , 20 , 0.0 ,  lc}; 
Point(8) = { 20 , 20 , 0.0 ,   lc}; 
Point(9) = { 20 , -20 , 0.0 ,  lc}; 

Line(5)={6,7};
Line(6)={7,8};
Line(7)={8,9};
Line(8)={9,6};

Line Loop(1)={5,6,7,8};
Line Loop(2)={1,2,3,4};

Plane Surface(1) = {2,1};

Coherence;

Field[1] = BoundaryLayer;
Field[1].EdgesList = {1,2,3,4};
Field[1].NodesList = {2,3,4,5};
Field[1].hfar = lc;
Field[1].hwall_n = 0.01;
Field[1].hwall_t = ht;
Field[1].thickness = 0.2;
Field[1].ratio = 1.1;
Field[1].Quads = 1;
Field[1].IntersectMetrics = 0;
BoundaryLayer Field = 1;


// Subdivision algorithm is not working properly for high-order boundary layer
Mesh.Algorithm = 9; // (0) Automatic (1) MeshAdapt (5) Delauney (6) Frontal (8) Delauney Quad (9) Packing of parallelograms
Mesh.RecombineAll=1;
Mesh.RecombinationAlgorithm = 0;  // (0) standard (1) blossom
Mesh.SubdivisionAlgorithm = 1; // (0) None (1) All Quads (2) All Hexas
Mesh.HighOrderOptimize = 1;

// Other Meshing options
Mesh.ElementOrder = 2;
Mesh.NumSubEdges = 10; // Show high-order elements smoothly
Mesh.RemeshAlgorithm = 0; // (0) no split (1) automatic (2) automatic only with metis
Mesh.RemeshParametrization = 0; // (0) harmonic (1) conformal spectral (7) conformal finite element
Mesh.Smoothing=0;

Mesh.HighOrderNumLayers = 10;
Mesh.HighOrderPoissonRatio = 0.33;
Mesh.HighOrderThresholdMin = 0.1;
Mesh.HighOrderThresholdMax = 2;
Mesh.HighOrderOptPrimSurfMesh = 0;

