/* ---------------------------------------------------------------------------------------------
This is a 3D finite element calculation for a large 1/4" thick steel plate
with small semi-spherical anomaly in the middle.
The plate is magnetized with two far strong magnets. A constant zero potential
surface is defined at one side of the anomaly an another at the other side (both
have the same zero potential). There is an air volume at each side of the plate.

Symetry is exploited to calculate just one quadrant of the problem.

Janvrot IVM,  Sep 2002
-----------------------------------------------------------------------------------------------*/

// Constants
// mm = 1e-3;

Include "par.inc";  // defect parameters

// Auxiliary values
sc = Sin(Pi/4);
re = (p * p + r * r) / (2 * p); // raio da superficie esferica

// Points for the boundary

Point( 2) = {0, -t, 0, clfd};
Point( 3) = {0, Rint, 0, clbp};
Point( 4) = {0, -Rint, 0, clbp};
Point(10) = {Rint, 0, 0, clbp};
Point(11) = {Rext, 0, 0, clbp};
Point(23) = {0, Rext, 0, clbp};
Point(24) = {0, -Rext, 0, clbp};


Point( 13) = {AcoX, 0, 0, clim};
Point( 14) = {AcoX, -t, 0, clim};

// Points for the anomaly

Point(1) = {0, 0, 0, clfd};
Point(17) = {r, 0, 0, clfd};
Point(20) = {0, -p, 0, clfd};
Point (21) = {0, re-p, 0, clfd}; // centro da elipse

// Points for the magnet

Point(33) = {AcoX + ImaX, 0, 0, clim};
Point(34) = {AcoX + ImaX, -t, 0, clim};

// Sufaces and volumes created with gmsh

Circle(4) = {20,21,17};
Line(5) = {4,2};
Line(6) = {2,20};
Line(7) = {20,3};
Line(10) = {33,34};
Line(13) = {34,14};
Line(14) = {14,2};
Line(15) = {33,13};
Line(16) = {13,17};
Line(17) = {13,14};


Circle(18) = {23,1,11};
Circle(19) = {11,1,24};
Circle(20) = {3,1,10};
Circle(21) = {10,1,4};
Line(22) = {23,3};
Line(23) = {10,11};
Line(24) = {4,24};
Line Loop(25) = {23,-18,22,20};
Plane Surface(26) = {25};
Line Loop(27) = {24,-19,-23,21};
Plane Surface(28) = {27};
Line Loop(29) = {5,-14,-13,-10,15,16,-4,7,20,21};
Plane Surface(30) = {29};
Line Loop(31) = {17,14,6,4,-16};
Plane Surface(32) = {31};
Line Loop(33) = {17,-13,-10,15};
Plane Surface(34) = {33};


// Boundary Physical Lines
Physical Line(LineA_Code) = {22,7,5,24};
Physical Line(LineC_Code) = {18,19};

// Physical Surface
Physical Surface(Steel_Code) = {32};
Physical Surface(Magnet_Code) = {34};
Physical Surface(Air_Code) = {30};
Physical Surface(AirInf_Code) = {26,28};
