% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.5a (COMSOL 3.5.0.608, $Date: 2009/05/11 07:38:49 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.5';
vrsn.ext = 'a';
vrsn.major = 0;
vrsn.build = 608;
vrsn.rcs = '$Name: v35ap $';
vrsn.date = '$Date: 2009/05/11 07:38:49 $';
fem.version = vrsn;

% Geometry
g1=cylinder3('0.4','0.1','pos',{'0','0','-0.05'},'axis',{'0','0','1'},'rot','0');

% Analyzed geometry
clear s
s.objs={g1};
s.name={'tissue'};
s.tags={'g1'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);
g2=cylinder3('0.0004','0.02','pos',{'0','0','0.03'},'axis',{'0','0','1'},'rot','0');

% Analyzed geometry
clear s
s.objs={g1,g2};
s.name={'tissue','pos'};
s.tags={'g1','g2'};

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

%new Alex negative electrode creation
E_number=4; %number of negative electrodes
R=0.02; % radius of a treated area
alpha=[]; % array of angels angel of electrode coordinate
a=0; % starting angel
for ii=1:E_number % creation of array of angels for negative electrodes
   alpha(ii)=a;
   a=a+2*pi/E_number;
end
for ii=1:E_number % x and y are arrays with electrode position coorinate In Comsol axis base point x and y. 
    x(ii)=R*cos(alpha(ii));
    y(ii)=R*sin(alpha(ii));
end
for ii=1:E_number
   gElectrode=cylinder3('0.0004','0.02','pos',{num2str(x(ii)),num2str(y(ii)),'0.03'},'axis',{'0','0','1'},'rot','0');
   s.objs{end+1}=gElectrode;
   s.name{end+1}=['negative'];
   %s.name{end+1}=['El' num2str(ii)];
   s.tags{end+1}=['El' num2str(ii)];
end
% end Alex new

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

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

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% (Default values are not included)

% end Alex new

% Application mode 1 for 
clear appl
appl.mode.class = 'EmConductiveMediaDC';
appl.module = 'ACDC';
appl.sshape = 2;
appl.border = 'on';
appl.assignsuffix = '_emdc';
clear bnd
bnd.V0 = {0,0,3000};

bnd.type = {'nJ0','V0','V'};

faces = domain2faces(fem);
 
 tissueDomain = 0;
 maxfaces=0;
 for i = 2 : length(faces)
    if (length(faces{i})>maxfaces)
        tissueDomain = i;
        maxfaces = length(faces{i});
    end;
end;
bnd.ind = [];
for i = 1 : length(faces)
    if (i==2)
        continue;
    end;
    
    for j = 1 : length(faces{i})
       
        bnd.ind(faces{i}(j)) = max(i-1,1);
    end;
end;


bnd.ind



appl.bnd = bnd;
clear equ
equ.sigma = {0.2,'10^7'};
equ.ind = [1,2];
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',{'V'}, ...
                  'outcomp',{'V'}, ...
                  'blocksize','auto', ...
                  'linsolver','cg', ...
                  'prefun','amg');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'slicedata',{'V','cont','internal','unit','V'}, ...
         'slicexspacing',5, ...
         'sliceyspacing',0, ...
         'slicezspacing',0, ...
         'slicemap','Rainbow', ...
         'title','Slice: Electric potential [V]', ...
         'grid','on', ...
         'campos',[-263.60143118283463,-343.532073434725,249.99999999999994], ...
         'camtarget',[0,0,0], ...
         'camup',[0,0,1], ...
         'camva',5.153143660537722);