0001 function [mapping, outside] = mk_approx_c2f( f_mdl, c_mdl );
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050 if ischar(f_mdl) && strcmp(f_mdl, 'UNIT_TEST'); do_unit_test; return; end
0051
0052 if ischar(f_mdl) && strcmp(f_mdl, 'LOAD'); load; return; end
0053
0054 [c_mdl, f_mdl] = assign_defaults( c_mdl, f_mdl );
0055
0056 copt.cache_obj = cache_obj(c_mdl, f_mdl);
0057 copt.fstr = 'mk_coarse_fine_mapping';
0058
0059
0060 [mapping, outside] = eidors_cache(@mapping_calc,{f_mdl,c_mdl},copt);
0061
0062
0063
0064 function [mapping, outside] = mapping_calc(f_mdl, c_mdl)
0065 f_mdl= offset_and_project( f_mdl, c_mdl);
0066 z_depth = c_mdl.mk_coarse_fine_mapping.z_depth;
0067
0068 f_elems = all_contained_elems( f_mdl, c_mdl, z_depth);
0069 mapping = contained_elems_i( f_mdl, c_mdl, f_elems, z_depth);
0070
0071 if isfield(c_mdl,'coarse2fine')
0072 mapping = mapping*c_mdl.coarse2fine;
0073 end
0074
0075
0076 if nargout>1;
0077 outside = 1 - sum(mapping,2);
0078 end
0079
0080
0081 function c_obj = cache_obj(c_mdl, f_mdl)
0082 c_obj = {c_mdl.nodes, c_mdl.elems, c_mdl.mk_coarse_fine_mapping, ...
0083 f_mdl.nodes, f_mdl.elems, f_mdl.interp_mesh};
0084
0085
0086
0087 function c_elems = all_contained_elems( fm, cm, z_depth)
0088 [nf,ef]= size(fm.elems);
0089 [nc,ec]= size(cm.elems);
0090 fm_pts = zeros(nf*ef,3);
0091
0092
0093
0094
0095 s_fac= .9;
0096 for dim= 1:ef-1
0097
0098 fm_pt= reshape(fm.nodes(fm.elems,dim),nf,ef);
0099 fm_ctr= mean(fm_pt,2)*ones(1,ef);
0100 fm_pt = s_fac*(fm_pt-fm_ctr) + fm_ctr;
0101 fm_pts(:,dim) = fm_pt(:);
0102 end
0103
0104 tsn= search_fm_pts_in_cm(cm, fm_pts, z_depth);
0105 tsn= reshape( tsn, [], ef);
0106
0107
0108
0109 c_elems= all(diff(tsn,1,2)==0,2) .* tsn(:,1);
0110 c_elems(all(tsn==-1,2))= -1;
0111
0112
0113
0114
0115
0116 function tsn= search_fm_pts_in_cm(cm, fm_pts, z_depth);
0117 dc= size(cm.elems,2)-1;
0118 df= size(fm_pts,2);
0119
0120 tsn= -ones(size(fm_pts,1),1);
0121 not_oor= (tsn==-1);
0122
0123 if dc==2
0124
0125 if df==3
0126
0127 not_oor= not_oor & any( abs(fm_pts(:,3) ) <= z_depth , 2);
0128 end
0129 dims=1:2;
0130
0131 elseif dc==3
0132
0133 dims=1:3;
0134
0135 else
0136 error('coarse model must be 2 or 3D');
0137 end
0138
0139 tsn(not_oor)= tsearchn(cm.nodes(:,dims), cm.elems, fm_pts(not_oor,dims));
0140 tsn(isnan(tsn))= -1;
0141
0142
0143
0144
0145
0146
0147 function c_elems = contained_elems_i( fm, cm, idx, z_depth)
0148 [nc,dc]= size(cm.elems);
0149 [nf,df]= size(fm.elems);
0150
0151 fidx= find(idx==0);
0152 l_fidx= length(fidx);
0153
0154 c_e_i= []; c_e_j=[]; c_e_v=[];
0155
0156
0157
0158 n_interp = 7-df;
0159 interp_mdl.nodes= fm.nodes;
0160 interp_mdl.elems= fm.elems(fidx,:);
0161 interp_mdl.interp_mesh.n_interp = fm.interp_mesh.n_interp;
0162
0163
0164 fm_pts = interp_mesh( interp_mdl, n_interp);
0165 l_interp = size(fm_pts,3);
0166
0167 fm_pts = permute( fm_pts, [3,1,2]);
0168 fm_pts = reshape(fm_pts, [], df-1);
0169
0170 tsn= search_fm_pts_in_cm(cm, fm_pts, z_depth);
0171
0172 tsn_idx= ones(l_interp,1)*fidx(:)';
0173 tsn_idx= tsn_idx(:);
0174
0175 outside_idx= tsn==-1;
0176 tsn(outside_idx) = [];
0177 tsn_idx(outside_idx) = [];
0178
0179 in_idx= find(idx<=0);
0180 ridx= 1:nf; ridx(in_idx)=[];
0181 idx(in_idx)=[];
0182
0183
0184
0185 c_elems = sparse(ridx,idx,1,nf,nc) + ...
0186 sparse(tsn_idx,tsn,1,nf,nc)/l_interp;
0187
0188
0189
0190 function xyz = interpxyz( xyzmin, xyzmax, n_interp)
0191 xyzdelta= xyzmax - xyzmin;
0192 xyz_interp = 1 + floor(n_interp * xyzdelta / max(xyzdelta) );
0193 xspace = linspace(xyzmin(1), xyzmax(1), xyz_interp(1) );
0194 yspace = linspace(xyzmin(2), xyzmax(2), xyz_interp(2) );
0195 zspace = linspace(xyzmin(3), xyzmax(3), xyz_interp(3) );
0196 [xx3,yy3,zz3] = ndgrid( xspace, yspace, zspace );
0197 xyz= [xx3(:), yy3(:), zz3(:)];
0198
0199
0200 function f_mdl= offset_and_project( f_mdl, c_mdl)
0201 [fn,fd]= size(f_mdl.nodes);
0202 T= c_mdl.mk_coarse_fine_mapping.f2c_offset;
0203 M= c_mdl.mk_coarse_fine_mapping.f2c_project;
0204
0205 f_mdl.nodes= (M*( f_mdl.nodes - ones(fn,1)*T )')';
0206
0207 function [c_mdl f_mdl] = assign_defaults( c_mdl, f_mdl )
0208 [fn,fd]= size(f_mdl.nodes);
0209 try
0210 c_mdl.mk_coarse_fine_mapping.f2c_offset;
0211 catch
0212 c_mdl.mk_coarse_fine_mapping.f2c_offset= zeros(1,fd);
0213 end
0214 try
0215 c_mdl.mk_coarse_fine_mapping.f2c_project;
0216 catch
0217 c_mdl.mk_coarse_fine_mapping.f2c_project= speye(fd);
0218 end
0219 try
0220 c_mdl.mk_coarse_fine_mapping.z_depth;
0221 catch
0222 c_mdl.mk_coarse_fine_mapping.z_depth= inf;
0223 end
0224 try
0225 f_mdl.interp_mesh.n_interp;
0226
0227
0228 catch
0229 f_mdl.interp_mesh.n_interp = 7 - size(f_mdl.elems,2);
0230 end
0231
0232
0233 function do_unit_test
0234 fmdl = mk_circ_tank(2,[],2); fmdl.nodes = fmdl.nodes*2;
0235 cmdl = mk_circ_tank(2,[],2); cmdl.nodes = cmdl.nodes*2;
0236 c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0237 unit_test_cmp('t1',c2f,eye(16),1e-15)
0238
0239 fmdl = mk_circ_tank(3,[],2);
0240 fmdl.nodes = fmdl.nodes*3;
0241 c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0242 unit_test_cmp('t2 analytic',c2f,[eye(16);zeros(20,16)],1e-15)
0243 c2f = mk_approx_c2f( fmdl, cmdl);
0244 unit_test_cmp('t2 approx',c2f,[eye(16);zeros(20,16)],1e-15)
0245
0246 fmdl = mk_circ_tank(2,[],2); fmdl.nodes = fmdl.nodes*2;
0247 cmdl = mk_circ_tank(1,[],2); cmdl.nodes = cmdl.nodes*1;
0248 c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0249 unit_test_cmp('t3 analytic',c2f,[eye(4);zeros(12,4)],1e-15)
0250 c2f = mk_approx_c2f( fmdl, cmdl);
0251 unit_test_cmp('t3 approx',c2f,[eye(4);zeros(12,4)],1e-15)
0252
0253 fac = 0.8;
0254 cmdl = mk_circ_tank(1,[],2); cmdl.nodes = cmdl.nodes*fac;
0255 c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0256 unit_test_cmp('t4 analytic',c2f,[eye(4)*fac^2;zeros(12,4)],1e-15)
0257 c2f = mk_approx_c2f( fmdl, cmdl);
0258 unit_test_cmp('t4 approx',c2f,[eye(4)*2/3;zeros(12,4)],1e-15)
0259
0260 fac=1.2;
0261 cmdl = mk_circ_tank(1,[],2); cmdl.nodes = cmdl.nodes*fac;
0262 c2f = mk_approx_c2f( fmdl, cmdl);
0263 unit_test_cmp('t5',c2f,[eye(4);eye(4)/3;kron(eye(4),[1;1])/15],1e-15);
0264
0265 fmdl = mk_circ_tank(10,[],2);
0266 cmdl = mk_circ_tank(8,[],2);
0267 c2f = mk_approx_c2f( fmdl, cmdl);
0268 unit_test_cmp('t6',sum(c2f'),ones(1,size(c2f,1)),1e-14);
0269
0270 cmdl.nodes = cmdl.nodes*0.95;
0271
0272 c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0273
0274 function load
0275
0276
0277 electrodes_per_plane = 16;
0278 number_of_planes = 2;
0279 tank_radius = 0.2;
0280 tank_height = 0.5;
0281 fine_mdl = ng_mk_cyl_models([tank_height,tank_radius],...
0282 [electrodes_per_plane,0.15,0.35],[0.01]);
0283
0284
0285 coarse_mdl_maxh = 0.07;
0286 coarse_mdl = ng_mk_cyl_models([tank_height,tank_radius,coarse_mdl_maxh],[0],[]);
0287
0288 disp('Calculating coarse2fine mapping ...');
0289 inv3d.fwd_model.coarse2fine = ...
0290 mk_coarse_fine_mapping( fine_mdl, coarse_mdl);
0291 disp(' ... done');