% 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  % create an array with negative electrodes
   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';


%**********Alex New 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=cell(size(domain));
domainNames{1}='surround';
domainNames{2}='tissue';
for i=3:length(domain)
    domainNames{i}='neg';
end

% find wich electrodeis positive
% if E_number is even then the positive electrode is the middle of
% electrode domains names 
if mod(E_number,2)==0
    domainNames{(3+length(domain))/2}='pos';
else
    % E_number is odd
    if mod(E_number+1,4)==0
        % if E_number is odd and E_number+1 is divisible by 4 then we should round up
        domainNames{ceil(E_number/2)+3}='pos';
    else
        % else we should round down
        domainNames{floor(E_number/2)+3}='pos';
    end
end
% define boundary conditions
clear bnd
bnd.V0 = {0,0,3000};
bnd.type = {'nJ0','V0','V'};

%faces = domain2faces(fem);

bnd.ind = ones(1,length(upDown));    % init all faces to be insulation. '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;    % 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}, 'pos'))
        bnd.ind(domain{i})=3;    %  initial condition for all faces on the positive electrode to be 3000V
    end
end



% 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;
% define permitivity
clear equ
equ.sigma = {0.2,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

%**********Alex New code                  
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);