// NAFEMS LE10 benchmark geometry & mesh
SetFactory("OpenCASCADE");

a = 1000;   // geometric parameters (in mm)
b = 2750;
c = 3250;
d = 2000;
h = 600;

Point(1) = {0, a, 0};  // define the four points A, B, C and D
Point(2) = {0, b, 0};
Point(3) = {c, 0, 0};
Point(4) = {d, 0, 0};


Line(1) = {1, 2};      // join them with the ellipses and straight edges
Ellipse (2) = {0,0,0, c, b, 0, Pi/2};
Line(3) = {3, 4};
Ellipse (4) = {0,0,0, d, a, 0, Pi/2};
Coherence;             // merge the points 

Curve Loop(1) = {1, -2, 3, 4};   // create the bottom surface
Plane Surface(1) = {1};

Extrude {0, 0, 0.5*h} { Surface{1}; }   // extrude it (twice to get the mid-plane edge)
Extrude {0, 0, 0.5*h} { Surface{6}; }

Physical Surface("DCD'C'") = {9, 4};    // define physical groups for boundary conditions
Physical Surface("ABA'B'") = {7, 2};
Physical Surface("BCB'C'") = {8, 3};
Physical Surface("upper") = {11};
Physical Curve("midplane") = {9};
Physical Volume("bulk") = {1,2};

// meshing settings
Mesh.ElementOrder = 2;          // use second-order 
Mesh.RecombineAll = 1;          // use hexas instead of tets
Mesh.SecondOrderLinear = 0;     // use curved elemetns
Mesh.SecondOrderIncomplete = 1; // use hex20 instead of hex27

// ask for a structured mesh
// the original 6x4x2 NAFEMS "fine" mesh is obtained with CharacteristicLengthFactor = 1
Transfinite Curve{2,9,17,4,12,20} = 6/Mesh.CharacteristicLengthFactor + 1;
Transfinite Curve{1,7,15,3,11,19} = 4/Mesh.CharacteristicLengthFactor + 1;
Transfinite Curve{5,13,6,14,8,16,10,18} = 1/Mesh.CharacteristicLengthFactor + 1;
Transfinite Surface {1,2,3,4,5,6,7,8,9,10,11};
Transfinite Volume {1,2};
