% 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.3','0.1','pos',{'0','0','-0.05'},'axis',{'0','0','1'},'rot','0');%Tissue cylinder

% 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');% Positive cylinder in the middle

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

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

N=4; % Number of negative electrodes
R=0.02; % Electrode array radius
alpha=0; % angel of electrode coordinate

for ii=1:N % creation of electrodes with negative potential
   gElectrode=cylinder3('0.0004','0.02','pos',{R*cos(alpha),R*sin(alpha),'0.03'},'axis',{'0','0','1'},'rot','0');
   alpha=alpha+2*pi/N;
   s.objs{end+1}=gElectrode;
   s.name{end+1}=['El' num2str(ii)];
   s.tags{end+1}=['El' num2str(ii)];
end

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

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

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','longest');

% (Default values are not included)

% Application mode 1
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','V','V'};
bnd.ind = [1,1,1,1,2,2,2,1,2,2,3,3,3,1,1,3,3,1];
appl.bnd = bnd;
clear equ
equ.sigma = {0.2,'5.5*10^7'};
equ.ind = [1,2,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',{'normE_emdc','cont','internal','unit','V/cm'}, ...
         'slicedlim',[500 80000], ...
         'slicexspacing',0, ...
         'sliceyspacing',0, ...
         'slicezspacing',5, ...
         'slicemap','Rainbow', ...
         'title','Slice: Electric field, norm [V/cm]', ...
         'grid','on', ...
         'campos',[-0.7798263573677466,-2.0932661536335715,4.473265883267742], ...
         'camtarget',[0,0,0], ...
         'camup',[0.4803155200573164,0.7592985067978727,0.43904758143803785], ...
         'camva',5.153143660537646);

% Integrate
I1=postint(fem,'normE_emdc>500[V/m]', ...
           'unit','m^3', ...
           'recover','off', ...
           'dl',1);

% Integrate
I2=postint(fem,'normE_emdc>80000[V/m]', ...
           'unit','m^3', ...
           'recover','off', ...
           'dl',1);