bld_master_full

function [Ef,D,Ela] = bld_master_full(vtx,simp,mat,elec,zc);

System matrix assembling based on the complete electrode model. 
This function is called within fem_master_full.



Ef   = The UNreferenced system matrix.
D    = The sgradients of the shape functions over each element.
Ela  = Normalised volums of the elements
vtx  = The vertices matrix. The coordinates of the nodes in 3D.
simp = The simplices matrix. Unstructured tetrahedral.
mat  = As for MATerial information. The conductivity vector.(isotropic)
elec = The matrix that holds the boundary faces assigned as electrodes. Its typical
       dimensions are [number of electrodes x 3*number of faces per electrode].
zc   = The array of electrode contact impedances.

0001 function [Ef,D,Ela] = bld_master_full(vtx,simp,mat,elec,zc);
0002 %function [Ef,D,Ela] = bld_master_full(vtx,simp,mat,elec,zc);
0003 %
0004 %System matrix assembling based on the complete electrode model.
0005 %This function is called within fem_master_full.
0006 %
0007 %
0008 %
0009 %Ef   = The UNreferenced system matrix.
0010 %D    = The sgradients of the shape functions over each element.
0011 %Ela  = Normalised volums of the elements
0012 %vtx  = The vertices matrix. The coordinates of the nodes in 3D.
0013 %simp = The simplices matrix. Unstructured tetrahedral.
0014 %mat  = As for MATerial information. The conductivity vector.(isotropic)
0015 %elec = The matrix that holds the boundary faces assigned as electrodes. Its typical
0016 %       dimensions are [number of electrodes x 3*number of faces per electrode].
0017 %zc   = The array of electrode contact impedances.
0018 
0019 dimen= size(vtx,2);
0020 if dimen==2
0021    [Ef,D,Ela] = bld_master_full_2d(vtx,simp,mat,elec,zc);
0022 elseif dimen==3
0023    [Ef,D,Ela] = bld_master_full_3d(vtx,simp,mat,elec,zc);
0024 else
0025    error('not 2d or 3d');
0026 end
0027 
0028 function [Ef,D,Ela] = bld_master_full_2d(vtx,simp,mat,elec,zc);
0029 
0030 [vr, vc] = size(vtx);
0031 [sr, sc] = size(simp);
0032 [er, ec] = size(elec);
0033 
0034 
0035 if length(mat) ~= sr
0036    error('Invalid conductivity information for this mesh');
0037 end
0038 
0039 
0040 if length(zc) == er
0041 
0042 
0043 %The column vector zc with the contact impedances in [Ohms] is required
0044 
0045 [E,D,Ela] = bld_master(vtx,simp,mat);
0046 
0047 
0048 E = full(E);
0049 
0050 Ef = spalloc(vr+er,vr+er, er * vr);
0051 
0052 Ef(1:vr,1:vr) = E;
0053 
0054 
0055 while length(zc) ~= er
0056       disp(sprintf('Please enter the correct zc column vector with length: %d',er));
0057       %[zc] = contact_impedance;
0058 end
0059 
0060 
0061 for q=1:er
0062    
0063    tang_dist = 0;
0064    
0065    q_th_ele = elec(q,:);  % Select the row of nodes corresponding to the current electrode
0066    
0067    q_th_ele_zf = nonzeros(q_th_ele)'; % Extract the dummy "zero" nodal numbers
0068    
0069    for w=1:2:length(q_th_ele_zf)
0070       
0071       m = q_th_ele_zf(w);
0072       n = q_th_ele_zf(w+1);
0073       
0074       % This way m & n nodes belong to the edge tangent to the electrode and also at the same simplex.
0075       
0076       % We now evaluate the distance "tangential contact" between m & n
0077       
0078       xm = vtx(m,1);
0079       ym = vtx(m,2); % m node coords
0080       xn = vtx(n,1);  
0081       yn = vtx(n,2); % n node coords
0082       
0083       [dist] = db2p(xm,ym,xn,yn); % distance mn
0084       
0085       cali_dist = dist ./ zc(q);  % coeficient for the distance mn
0086       
0087       tang_dist = tang_dist + cali_dist;
0088       
0089       % Start modifying "expanding" the E master matrix
0090       
0091       Ef(m,vr+q) = Ef(m,vr+q) - cali_dist/2 ; % Kv -> Ec  -> Vertical bar
0092       Ef(n,vr+q) = Ef(n,vr+q) - cali_dist/2 ; % Kv -> Ec
0093       
0094       Ef(vr+q,m) = Ef(vr+q,m) - cali_dist/2 ; % Kv' -> Ec' -> Horizontal bar
0095       Ef(vr+q,n) = Ef(vr+q,n) - cali_dist/2 ; % Kv' -> Ec'
0096       
0097       
0098       Ef(m,m) = Ef(m,m) + cali_dist/3; % Kz -> E -> Main bar
0099       Ef(n,n) = Ef(n,n) + cali_dist/3; % Kz -> E
0100       Ef(m,n) = Ef(m,n) + cali_dist/6; % Kz -> E
0101       Ef(n,m) = Ef(n,m) + cali_dist/6; % Kz -> E
0102       
0103    end % dealing with this electrode
0104    
0105    Ef(vr+q,vr+q) = Ef(vr+q,vr+q) + tang_dist;
0106    
0107 end %for the whole set of electrodes
0108 
0109 end
0110 
0111 function [Ef,D,Ela] = bld_master_full_3d(vtx,simp,mat,elec,zc);
0112 [vr, vc] = size(vtx);
0113 [sr, sc] = size(simp);
0114 [er, ec] = size(elec);
0115 
0116 
0117 if length(mat) ~= sr
0118    error('Invalid conductivity information for this mesh');
0119 end
0120 
0121 
0122 [Ef,D,Ela] = bld_master(vtx,simp,mat);
0123 
0124 
0125 % Add zeros so Ef is of size (vr+er) x (vr+er)
0126 [Ef_i, Ef_j, Ef_s]= find( Ef );
0127 Ef = sparse(Ef_i, Ef_j, Ef_s, vr+er, vr+er);
0128 
0129 
0130 %Up to this point we have calculated the master matrix without the influence of contact impedance.
0131 
0132 %The column vector zc with the contact
0133 %impedances in [Ohms] is required
0134 if length(zc) ~= er
0135       error(sprintf('zc (=%d) should be equal to er (=%d)',length(zc),er));
0136 end
0137 
0138 
0139 for q=1:er
0140    
0141    tang_area = 0;
0142    
0143    q_th_ele = nonzeros(elec(q,:));  % Select the row of nodes corresponding to the current electrode
0144    
0145    if length(q_th_ele) ==1 % check if point electrode
0146       m = q_th_ele;
0147       cali_area = 1 / zc(q);
0148    
0149       tang_area = tang_area + cali_area;
0150       
0151       Ef(m,vr+q) = Ef(m,vr+q) - cali_area/2 ; 
0152       Ef(vr+q,m) = Ef(vr+q,m) - cali_area/2 ; 
0153       
0154       Ef(m,m) = Ef(m,m) + cali_area/2;
0155 
0156    else % not point electrode - use complete electrode model
0157    for w=1:3:length(q_th_ele)
0158       
0159       m = q_th_ele(w);
0160       n = q_th_ele(w+1);
0161       l = q_th_ele(w+2);
0162       
0163         
0164       % This way m & n nodes belong to the edge tangential to the electrode
0165       % and also at the same simplex.
0166       % We will now evaluate the distance "tangential contact area" between m,n & l
0167       Are = triarea3d(vtx([m n l],:));
0168           
0169     % area mnl
0170       
0171       cali_area = (2*Are) ./ zc(q);  % coefficient for the area mnl
0172       %|J_k| = 2*Are
0173       
0174       tang_area = tang_area + cali_area;
0175       
0176       % Start modifying "expanding" the E master matrix
0177       
0178       Ef(m,vr+q) = Ef(m,vr+q) - cali_area/6 ; % Kv -> Ec  -> Vertical bar
0179       Ef(n,vr+q) = Ef(n,vr+q) - cali_area/6 ; 
0180       Ef(l,vr+q) = Ef(l,vr+q) - cali_area/6 ;
0181             
0182       Ef(vr+q,m) = Ef(vr+q,m) - cali_area/6 ; % Kv' -> Ec' -> Horizontal bar
0183       Ef(vr+q,n) = Ef(vr+q,n) - cali_area/6 ; 
0184       Ef(vr+q,l) = Ef(vr+q,l) - cali_area/6 ;
0185       
0186       Ef(m,m) = Ef(m,m) + cali_area/12; % Kz -> E -> Main bar
0187       Ef(m,n) = Ef(m,n) + cali_area/24;       
0188       Ef(m,l) = Ef(m,l) + cali_area/24;
0189       
0190       Ef(n,m) = Ef(n,m) + cali_area/24;
0191       Ef(n,n) = Ef(n,n) + cali_area/12; 
0192       Ef(n,l) = Ef(n,l) + cali_area/24;
0193       
0194       Ef(l,m) = Ef(l,m) + cali_area/24;
0195       Ef(l,n) = Ef(l,n) + cali_area/24;
0196       Ef(l,l) = Ef(l,l) + cali_area/12;
0197     
0198       
0199    end % dealing with this electrode
0200    end % point electrode
0201    Ef(vr+q,vr+q) = Ef(vr+q,vr+q) + 0.5*tang_area;
0202    
0203 end %for the whole set of electrodes
0204 
0205 % calculate distance between two points
0206 function [dist] = db2p(xa,ya,xb,yb);
0207 
0208    dist = sqrt((xb - xa).^2 + (yb - ya).^2);
0209 
0210 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0211 % This is part of the EIDORS suite.
0212 % Copyright (c) N. Polydorides 2003
0213 % Copying permitted under terms of GNU GPL
0214 % See enclosed file gpl.html for details.
0215 % EIDORS 3D version 2.0
0216 % MATLAB version 5.3 R11
0217 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%