% 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)


%*****************************************************************
% begin edited code
%*****************************************************************

% determine what faces are associated with which domains
% get the up/down info
upDown=geominfo(fem.geom,'out','ud');
% domain(n) holds the faces associated with domain n-1
% domain 0 is the volume surrounding the model
domain=cell(1,max(upDown(:))+1);
% loop through the faces and add them to the appropriate domains
for i=1:length(upDown)
    domain{upDown(1,i)+1}=[domain{upDown(1,i)+1},i];
    domain{upDown(2,i)+1}=[domain{upDown(2,i)+1},i];
end

% determine the names associated with the domains
domainNames={'surround','tissue','neg','neg','pos','neg','neg'};

% define boundary conditions
clear bnd
bnd.V0 = {0,0,3000};    % define negative voltage as 0V and positive voltage as 3000V
bnd.type = {'fp','V','V'};    % was complaining about type nJ0 so I changed it to a floating potential so it would compile
bnd.ind = ones(1,length(upDown));    % init all faces to be floating potential. 'length(upDown)' returns the number of faces

% find the faces that are negative
for i=1:length(domainNames)
    if (strcmp(domainNames{i}, 'neg'))
        bnd.ind(domain{i})=2;    % change the initial condition for all faces on all negative electrodes to be 0V
    end
end

% find the faces that are positive
for i=1:length(domainNames)
    if (strcmp(domainNames{i}, 'tip'))
        bnd.ind(domain{i})=3;    % change the initial condition for all faces on the positive electrode to be 3000V
    end
end

% define permitivity
clear equ
equ.sigma = {0.2,5.5e7};    % 5.5*10^7 is better written as 5.5e7
equ.ind = [];     % initialise the indices vector. when done, this vector will define which of the two
                  % options you have defined for each subdomain, sigma=0.2 or sigma=5.5e7

for i=1:length(domainNames)
    switch domainNames{i}
        case {'neg' 'pos'}
            equ.ind=[equ.ind,2];
        case {'tissue'}
            equ.ind=[equ.ind,1];
    end
end

% initialise the app1 variable
% Application mode 1
clear appl
appl.mode.class = 'EmConductiveMediaDC';
appl.module = 'ACDC';
appl.sshape = 2;
appl.border = 'on';
appl.assignsuffix = '_emdc';

% boundary conditions
appl.bnd = bnd;
% subdoiman conditions
appl.equ = equ;

%*****************************************************************
% end edited code
%*****************************************************************


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);