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