smp.geo

Include "smp_data.geo";

SetFactory("OpenCASCADE");

Mesh.Optimize = 1;

// circular x-section coil
vC()+=newv; Cylinder(newv) = {0,-Ly/2,0, 0,Ly,0, RintCoil};
vC()+=newv; Cylinder(newv) = {0,-Ly/2,0, 0,Ly,0, RextCoil};

vCoil()+=newv; BooleanDifference(newv) = { Volume{vC(1)}; Delete; }{ Volume{vC(0)}; Delete; };

// Smp 
vCoreE()=newv; Cylinder(newv) = { 0,-H/2,0,  0,H,0,  R};

// Air around
vA()+=newv; Sphere(newv) = {0,0,0, Rint};
vA()+=newv; Sphere(newv) = {0,0,0, Rext};

vAir()+=newv; BooleanDifference(newv) = { Volume{vA(1)}; Delete; }{ Volume{vA(0)}; };
vAir()+=newv; BooleanDifference(newv) = { Volume{vA(0)}; Delete; }{ Volume{vCoreE(),vCoil()}; };

If(Flag_Symmetry)
  xa =  -Rext;   ya =  -Rext;  za =  -Rext;
  dxa =  2*Rext; dya = 2*Rext; dza = 2*Rext;

  If(Flag_Symmetry==1)
    vAux=newv; Box(newv) = {xa, ya, 0*za, dxa, dya, -dza/2};
  EndIf
  If(Flag_Symmetry==2)
    vAux=newv; Box(newv) = {0*xa, ya, 0*za, dxa/2, dya, -dza/2};
  EndIf

  For k In {0:#vCoil()-1}
    vvv=newv; BooleanIntersection(newv) = { Volume{vCoil(k)}; Delete; }{ Volume{vAux}; };
    vCoil(k) = vvv;
  EndFor
  For k In {0:#vAir()-1}
    vvv=newv; BooleanIntersection(newv) = { Volume{vAir(k)}; Delete; }{ Volume{vAux}; };
    vAir(k) = vvv;
  EndFor

  vvv=newv; BooleanIntersection(newv) = { Volume{vCoreE(0)}; Delete; }{ Volume{vAux}; Delete; };
  vCoreE(0) = vvv;
EndIf

BooleanFragments{ Volume{vAir(), vCoreE(), vCoil()}; Delete; }{} // This needs to be done at the end

// Adapting mesh size
// characteristic lengths
lc0 = Pi*Rint/5; // air
lc1  = wcoil/3; // coil
lc2  = R/2; // smp

Characteristic Length { PointsOf{ Volume{vAir(0)}; } } = lc0;
Characteristic Length { PointsOf{ Volume{vCoil()}; } } = lc1;
Characteristic Length { PointsOf{ Volume{vCoreE()}; } } = lc2;


// Boundary conditions
tol = 1e-5;
cut_xy() = Surface In BoundingBox {-Rext-tol,-Rext-tol,-tol, 2*(Rext+tol), 2*(Rext+tol), 2*tol}; // 1/2, 1/4
cut_yz() = Surface In BoundingBox {-tol,-Rext-tol,-Rext-tol, 2*tol, 2*(Rext+tol), 2*(Rext+tol)}; // 1/4

all_vol() = Volume '*';
Printf("all_vol()=", all_vol());
bndAir() = CombinedBoundary{Volume{all_vol()};};
bndAir() -= {cut_xy(),cut_yz()};


//=================================================
// Some colors... for aesthetics :-)
//=================================================
Recursive Color SkyBlue {Volume{vAir()};}
Recursive Color SteelBlue {Volume{vCoreE()};}
Recursive Color Red {Volume{vCoil()};}

//=================================================
// Physical regions for FE analysis with GetDP
//=================================================

Physical Volume(ECORE) = vCoreE();
bnd_vCoreE() = CombinedBoundary{Volume{vCoreE()};};
bnd_vCoreE() -= {cut_xy(),cut_yz()};
Physical Surface(SKINECORE) = bnd_vCoreE();

Physical Surface(REFCORE)= bnd_vCoreE(1);
Physical Surface(REFCOREBOTTOM)= bnd_vCoreE(2);
pts() = PointsOf {Surface { bnd_vCoreE(1) }; };
Physical Point(REFCOREPOINT)= pts(0);

Physical Volume(COIL) = vCoil();
bnd_vCoil() = CombinedBoundary{Volume{vCoil()};};
bnd_vCoil() -= {cut_xy(),cut_yz()};
Physical Surface(SKINCOIL) = bnd_vCoil();

If(Flag_Infinity==0)
  Physical Volume(AIR) = {vAir()};
EndIf
If(Flag_Infinity==1)
  Physical Volume(AIR) = vAir(1);
  Physical Volume(AIRINF) = vAir(0);
EndIf
Physical Surface(SURF_AIROUT) = bndAir();

Physical Surface(CUT_XY) = {cut_xy()};
Physical Surface(CUT_YZ) = {cut_yz()}; // BC if symmetry