0001 function [Ef,D,Ela] = bld_master_full(vtx,simp,mat,elec,zc);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
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
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
0060 end
0061
0062
0063 for q=1:er
0064
0065 tang_dist = 0;
0066
0067 q_th_ele = elec(q,:);
0068
0069 q_th_ele_zf = nonzeros(q_th_ele)';
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
0077
0078
0079
0080 xm = vtx(m,1);
0081 ym = vtx(m,2);
0082 xn = vtx(n,1);
0083 yn = vtx(n,2);
0084
0085 [dist] = db2p(xm,ym,xn,yn);
0086
0087 cali_dist = dist ./ zc(q);
0088
0089 tang_dist = tang_dist + cali_dist;
0090
0091
0092
0093 Ef(m,vr+q) = Ef(m,vr+q) - cali_dist/2 ;
0094 Ef(n,vr+q) = Ef(n,vr+q) - cali_dist/2 ;
0095
0096 Ef(vr+q,m) = Ef(vr+q,m) - cali_dist/2 ;
0097 Ef(vr+q,n) = Ef(vr+q,n) - cali_dist/2 ;
0098
0099
0100 Ef(m,m) = Ef(m,m) + cali_dist/3;
0101 Ef(n,n) = Ef(n,n) + cali_dist/3;
0102 Ef(m,n) = Ef(m,n) + cali_dist/6;
0103 Ef(n,m) = Ef(n,m) + cali_dist/6;
0104
0105 end
0106
0107 Ef(vr+q,vr+q) = Ef(vr+q,vr+q) + tang_dist;
0108
0109 end
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
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
0133
0134
0135
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,:));
0146
0147 if length(q_th_ele) ==1
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
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
0167
0168
0169 Are = triarea3d(vtx([m n l],:));
0170
0171
0172
0173 cali_area = (2*Are) ./ zc(q);
0174
0175
0176 tang_area = tang_area + cali_area;
0177
0178
0179
0180 Ef(m,vr+q) = Ef(m,vr+q) - cali_area/6 ;
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 ;
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;
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
0202 end
0203 Ef(vr+q,vr+q) = Ef(vr+q,vr+q) + 0.5*tang_area;
0204
0205 end
0206
0207
0208 function [dist] = db2p(xa,ya,xb,yb);
0209
0210 dist = sqrt((xb - xa).^2 + (yb - ya).^2);
0211
0212
0213
0214
0215
0216
0217
0218
0219