% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.5a (COMSOL 3.5.0.603, $Date: 2008/12/03 17:02:19 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.5';
vrsn.ext = 'a';
vrsn.major = 0;
vrsn.build = 603;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2008/12/03 17:02:19 $';
fem.version = vrsn;

Inv=1/100;
num=80;
Detz=Inv/num;
z_i=-0.0185;
z_sus_vec=(z_i-(Inv/2):Detz:z_i+(Inv/2));

for i=1:num
    
flclear fem

z_sus_vec_new(i)=z_sus_vec(i);

% Geometry
g1=cylinder3('0.0381','0.0127','pos',{'0','0','-0.00635'},'axis',{'0','0','1'},'rot','0');
g2=cylinder3('0.01905','0.00635','pos',{'0','0',z_sus_vec(i)},'axis',{'0','0','1'},'rot','0');
g3=block3('0.15','0.15','0.15','base','center','pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0');

% Analyzed geometry
clear s
s.objs={g1,g2,g3};
s.name={'CYL1','CYL2','BLK1'};
s.tags={'g1','g2','g3'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',4, ...
                  'hmaxsub',[2,0.002,3,0.001]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',4, ...
                  'hmaxsub',[2,0.002,3,0.001]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'MagnetostaticsNoCurrents';
appl.module = 'ACDC';
appl.sshape = 2;
appl.assignsuffix = '_emnc';
clear bnd
bnd.type = {'nB0','cont'};
bnd.ind = [1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,1];
appl.bnd = bnd;
clear equ
equ.magconstrel = {'mur','M','M'};
equ.M = {{0;0;0},{0;0;186211},{0;0;166078}};
equ.maxwell = {{},{},'sus'};
equ.ind = [1,2,3];
appl.equ = equ;
fem.appl{1} = appl;
fem.frame = {'ref'};
fem.border = 1;
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;

% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femstatic(fem, ...
                  'solcomp',{'Vm'}, ...
                  'outcomp',{'Vm'}, ...
                  'blocksize','auto', ...
                  'linsolver','cg', ...
                  'prefun','amg');

% Save current fem structure for restart purposes
data=postglobaleval(fem,{'sus_forcez_emnc'});
z_force_vec(i)=data.y;
end
plot(z_sus_vec_new,z_force_vec);