Description of bld

0001 function [Ef,D,Ela] = bld_master(vtx,simp,mat_ref);
0002 %function [Ef,D,Ela] = bld_master(vtx,simp,mat_ref);
0003 %
0004 %Builds up the main compartment (GAP-SHUNT) of the system matrix
0005 %for the complete electrode model. It is called within the function
0006 %fem_master_full.
0007 %
0008 %
0009 %
0010 %Ef      = The UNreferenced GAP based system matrix
0011 %D       = The sgradients of the shape functions
0012 %          over each element.
0013 %Ela     = Normalised volums of the elements
0014 %vtx     = The vertices matrix.
0015 %simp    = The simplices matrix.
0016 %mat_ref = The reference CONDUCTIVITY at each element.
0017 %In the complex case mat_ref(i) = sigma(i) - epsilon(i)
0018 
0019 warning('EIDORS:deprecated','BLD_MASTER is deprecated as of 07-Jun-2012. ');
0020 
0021 dimen= size(vtx,2);
0022 if dimen==2
0023    [Ef,D,Ela] = bld_master_2d(vtx,simp,mat_ref); 
0024 elseif dimen==3
0025    [Ef,D,Ela] = bld_master_3d(vtx,simp,mat_ref);
0026 else
0027    error('not 2d or 3d');
0028 end
0029    
0030 function [Ef,A,Ela] = bld_master_2d(vtx,simp,mat_ref) 
0031 %Based on gm_assemble by S.Vavasis , http://www.cs.cornel.edu/
0032 
0033 [vr, vc] = size(vtx);
0034 [sr, sc] = size(simp);
0035 a = mat_ref;
0036 
0037 ilist = kron((1:vc*sr), [1,1]);
0038 jlist = zeros(1,sr*vc*2);
0039 slist = zeros(1,sr*vc*2);
0040 
0041 for d = 1 : vc
0042   jlist(2 * d - 1 : 2 * vc : sr * vc * 2) = simp(:,1)';
0043   jlist(2 * d : 2 * vc : sr * vc * 2) = simp(:, d + 1);
0044   slist(2 * d - 1 : 2 * vc : sr * vc * 2) = -ones(1,sr);
0045   slist(2 * d : 2 * vc : sr * vc * 2) = ones(1,sr);
0046 end
0047 
0048 A0 = sparse(ilist,jlist,slist,vc*sr,vr);
0049 
0050 
0051 if vc == 2
0052   J1 = A0 * vtx(:,1);
0053   J2 = A0 * vtx(:,2);
0054   J11 = J1(1:2:sr*2);
0055   J12 = J1(2:2:sr*2);
0056   J21 = J2(1:2:sr*2);
0057   J22 = J2(2:2:sr*2);
0058   detJ = J11 .* J22 - J21 .* J12;
0059   
0060   
0061   invJ11 = J22 ./ detJ;
0062   invJ12 = -J12 ./ detJ;
0063   invJ21 = -J21 ./ detJ;
0064   invJ22 = J11 ./ detJ;
0065 elseif vc == 3
0066   J1 = A0 * vtx(:,1);
0067   J2 = A0 * vtx(:,2);
0068   J3 = A0 * vtx(:,3);
0069   J11 = J1(1:3:sr*3);
0070   J12 = J1(2:3:sr*3);
0071   J13 = J1(3:3:sr*3);
0072   J21 = J2(1:3:sr*3);
0073   J22 = J2(2:3:sr*3);
0074   J23 = J2(3:3:sr*3);
0075   J31 = J3(1:3:sr*3);
0076   J32 = J3(2:3:sr*3);
0077   J33 = J3(3:3:sr*3);
0078   detJ = J11 .* J22 .* J33 - J11 .* J23 .* J32 - J12 .* J21 .* J33 ...
0079           + J12 .* J23 .* J31 + J13 .* J21 .* J32 - J13 .* J22 .* J31;
0080        
0081        
0082   invJ11 = (J22 .* J33 - J23 .* J32) ./ detJ;
0083   invJ12 = (J32 .* J13 - J12 .* J33) ./ detJ;
0084   invJ13 = (J12 .* J23 - J22 .* J13) ./ detJ;
0085   invJ21 = (J31 .* J23 - J21 .* J33) ./ detJ;
0086   invJ22 = (J11 .* J33 - J13 .* J31) ./ detJ;
0087   invJ23 = (J21 .* J13 - J11 .* J23) ./ detJ;
0088   invJ31 = (J21 .* J32 - J31 .* J22) ./ detJ;
0089   invJ32 = (J31 .* J12 - J11 .* J32) ./ detJ;
0090   invJ33 = (J11 .* J22 - J21 .* J12) ./ detJ;
0091 else
0092   error('Master matrix construction failed')
0093 end
0094 
0095 
0096 ilist = kron((1 : vc * sr), ones(1,vc));
0097 jlist = zeros(1, sr*vc^2);
0098 for d = 1 : vc 
0099   jlist(d:vc:sr*vc^2) = kron((d:vc:vc*sr),ones(1,vc));
0100 end
0101 
0102 if vc == 2
0103   slist = zeros(1,sr*4);
0104   slist(1:4:sr*4) = invJ11;
0105   slist(2:4:sr*4) = invJ21;
0106   slist(3:4:sr*4) = invJ12;
0107   slist(4:4:sr*4) = invJ22;
0108 else
0109   slist = zeros(1,sr*9);
0110   slist(1:9:sr*9) = invJ11;
0111   slist(2:9:sr*9) = invJ21;
0112   slist(3:9:sr*9) = invJ31;
0113   slist(4:9:sr*9) = invJ12;
0114   slist(5:9:sr*9) = invJ22;
0115   slist(6:9:sr*9) = invJ32;
0116   slist(7:9:sr*9) = invJ13;
0117   slist(8:9:sr*9) = invJ23;
0118   slist(9:9:sr*9) = invJ33;
0119 end
0120 
0121 
0122 ElJac = sparse(ilist,jlist,slist,vc*sr,vc*sr);
0123 A = ElJac * A0;
0124 
0125 
0126 Vols = abs(detJ) / prod(1:vc);
0127 
0128 materials = length(a);
0129 volumes = size(Vols);
0130 
0131 if materials ~= volumes
0132   error('Some elements have not been assigned');
0133 end;
0134 
0135 Ela = sparse( (1:vc*sr), (1:vc*sr),kron( (a .* Vols)',ones(1,vc)) );
0136 
0137 Ef = A'*Ela*A; 
0138 
0139 %This is for the Jacobian matrix (does not include conductivity)
0140 Ela = sparse( (1:vc*sr), (1:vc*sr),kron( Vols.',ones(1,vc)) );
0141 
0142 
0143 %
0144 % 3D BLD MASTER
0145 %
0146 function [Ef,D,Ela] = bld_master_3d(vtx,simp,mat_ref) 
0147 [nv, dimen] = size(vtx);  %Number of vertices and dimension
0148 [ns, dimen_p1] = size(simp); %Number of simplices
0149 a = mat_ref;
0150 dimen2 = 2*dimen;
0151 
0152 ils = 1:dimen*ns;
0153 ilst(1:2:dimen2*ns) = ils;
0154 ilst(2:2:dimen2*ns) = ils;
0155 
0156 
0157 patt = 1:dimen2:ns*dimen2;
0158 
0159 jlst(patt) = simp(:,1);
0160 jlst(patt+1) = simp(:,2);
0161 jlst(patt+2) = simp(:,1);
0162 jlst(patt+3) = simp(:,3);
0163 jlst(patt+4) = simp(:,1);
0164 jlst(patt+5) = simp(:,4);
0165 
0166 
0167 sls = ones(1,dimen*ns);
0168 slst(1:2:dimen2*ns) = -sls;
0169 slst(2:2:dimen2*ns) = sls;
0170 
0171 
0172 D0 = sparse(ilst,jlst,slst,dimen*ns,nv); 
0173 %D0 is the matrix of the definitions of the gradients on elements
0174 
0175 J1 = D0 * vtx(:,1);
0176 J2 = D0 * vtx(:,2);
0177 J3 = D0 * vtx(:,3);
0178   
0179 JJ= zeros( 3,3, ns );
0180 for w=1:3
0181    r=1;
0182    JJ(r,w,:)   = J1(w:dimen:ns*dimen);
0183    JJ(r+1,w,:) = J2(w:dimen:ns*dimen);
0184    JJ(r+2,w,:) = J3(w:dimen:ns*dimen);
0185 end
0186 
0187 dj = squeeze(sum( [prod([JJ(1,1,:);JJ(2,2,:);JJ(3,3,:)]); prod([JJ(1,2,:);JJ(2,3,:);JJ(3,1,:)]);...
0188                    prod([JJ(1,3,:);JJ(2,1,:);JJ(3,2,:)]); prod([-JJ(1,3,:);JJ(2,2,:);JJ(3,1,:)]);...
0189                    prod([-JJ(1,1,:);JJ(2,3,:);JJ(3,2,:)]); prod([-JJ(1,2,:);JJ(2,1,:);JJ(3,3,:)])]));
0190 
0191       
0192 ilst = kron((1:dimen*ns), ones(1,dimen));
0193 jlst = zeros(1, ns*dimen^2);
0194 for d = 1:dimen 
0195   jlst(d:dimen:ns*dimen^2) = kron((d:dimen:dimen*ns),ones(1,dimen));
0196 end
0197 slst = zeros(1,ns*dimen^2); 
0198    
0199   pat = 1:dimen^2:ns*dimen^2; 
0200   
0201   slst(pat) = squeeze(sum([prod([JJ(2,2,:) ; JJ(3,3,:)]); prod([-JJ(2,3,:) ; JJ(3,2,:)])])) ./ dj; 
0202   slst(pat+1) = squeeze(sum([prod([JJ(3,1,:) ; JJ(2,3,:)]); prod([-JJ(2,1,:) ; JJ(3,3,:)])])) ./ dj; 
0203   slst(pat+2) = squeeze(sum([prod([JJ(2,1,:) ; JJ(3,2,:)]); prod([-JJ(3,1,:) ; JJ(2,2,:)])])) ./ dj; 
0204   slst(pat+3) = squeeze(sum([prod([JJ(3,2,:) ; JJ(1,3,:)]); prod([-JJ(1,2,:) ; JJ(3,3,:)])])) ./ dj; 
0205   slst(pat+4) = squeeze(sum([prod([JJ(1,1,:) ; JJ(3,3,:)]); prod([-JJ(1,3,:) ; JJ(3,1,:)])])) ./ dj; 
0206   slst(pat+5) = squeeze(sum([prod([JJ(3,1,:) ; JJ(1,2,:)]); prod([-JJ(1,1,:) ; JJ(3,2,:)])])) ./ dj; 
0207   slst(pat+6) = squeeze(sum([prod([JJ(1,2,:) ; JJ(2,3,:)]); prod([-JJ(2,2,:) ; JJ(1,3,:)])])) ./ dj; 
0208   slst(pat+7) = squeeze(sum([prod([JJ(2,1,:) ; JJ(1,3,:)]); prod([-JJ(1,1,:) ; JJ(2,3,:)])])) ./ dj; 
0209   slst(pat+8) = squeeze(sum([prod([JJ(1,1,:) ; JJ(2,2,:)]); prod([-JJ(2,1,:) ; JJ(1,2,:)])])) ./ dj; 
0210   
0211 
0212 LocJac = sparse(ilst,jlst,slst,dimen*ns,dimen*ns);
0213 
0214 D = LocJac * D0;
0215 
0216 
0217 Vols = abs(dj(:)) / 6;
0218 
0219 materials = length(a);
0220 volumes = size(Vols);
0221 
0222 if materials ~= volumes
0223   error('Some elements have not been assigned');
0224 end;
0225 
0226 %This is for the global conductance matrix (includes conductivity)
0227 Ela = sparse( (1:dimen*ns), (1:dimen*ns),kron( (a.* Vols).',ones(1,dimen)) );
0228 
0229 Ef = D'*Ela*D; 
0230 
0231 %This is for the Jacobian matrix (does not include conductivity)
0232 Ela = sparse( (1:dimen*ns), (1:dimen*ns),kron( Vols.',ones(1,dimen)) );
0233 
0234 
0235 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0236 % This is part of the EIDORS suite.
0237 % Copyright (c) N. Polydorides 2003
0238 % Copying permitted under terms of GNU GPL
0239 % See enclosed file gpl.html for details.
0240 % EIDORS 3D version 2.0
0241 % MATLAB version 5.3 R11
0242 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0243
bld_master

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
