mk_approx_c2f

 MK_APPROX_C2F: create a mapping matrix from coarse to fine FEM
 [c2f,out]= mk_coarse_fine_mapping( f_mdl, c_mdl );
  
 Parameters:
    c_mdl is coarse fwd_model (no holes, see warning below)
    f_mdl is fine fwd_model

 C2F_ij is the fraction if f_mdl element i which is
   contained in c_mdl element j. This is used to map
   from data on the reconstruction model (c_mdl) to
   the forward model f_mdl as 
      elem_data_fine = Mapping*elem_data_coase

 OUT_i is the fraction of f_mdl element i which is not
   contained in any c_mdl element.

 OPTIONS:
 if the geometry of the fine and coarse models are not
  aligned, then they can be translated and mapped using
    coarse_xyz = (M* (fine_xyz - T)')'
  where
    T= c_mdl.mk_coarse_fine_mapping.f2c_offset (1xN_dims)
    M= c_mdl.mk_coarse_fine_mapping.f2c_project (N_dimsxN_dims)
  by default T= [0,0,0] and M=eye(3)

 if c_mdl is 2D and f_mdl is 3D, then parameter
     c_mdl.mk_coarse_fine_mapping.z_depth
     indicates the +/- z_depth which elements in 2D are
     considered to be extruded in 3D (default inf)

 NOTES:
 if c_mdl and f_mdl do not cover the same area, then 
    sum(c2f') will not be 1. If all coarse elements cover
    at least partially the fine ones, then this 
    can be corrected by:
      c2f = c2f./(sum(c2f,2) * ones(1,size(c2f,2)));
 if not all coarse elements cover fine ones, then this
    approach cannot be used. This will be fixed in a 
    future release

 WARNING:
 If c_mdl is not simply connected, the results are wrong!

 See also MK_C2F_CIRC_MAPPING

0001 function [mapping, outside] = mk_approx_c2f( f_mdl, c_mdl );
0002 % MK_APPROX_C2F: create a mapping matrix from coarse to fine FEM
0003 % [c2f,out]= mk_coarse_fine_mapping( f_mdl, c_mdl );
0004 %
0005 % Parameters:
0006 %    c_mdl is coarse fwd_model (no holes, see warning below)
0007 %    f_mdl is fine fwd_model
0008 %
0009 % C2F_ij is the fraction if f_mdl element i which is
0010 %   contained in c_mdl element j. This is used to map
0011 %   from data on the reconstruction model (c_mdl) to
0012 %   the forward model f_mdl as
0013 %      elem_data_fine = Mapping*elem_data_coase
0014 %
0015 % OUT_i is the fraction of f_mdl element i which is not
0016 %   contained in any c_mdl element.
0017 %
0018 % OPTIONS:
0019 % if the geometry of the fine and coarse models are not
0020 %  aligned, then they can be translated and mapped using
0021 %    coarse_xyz = (M* (fine_xyz - T)')'
0022 %  where
0023 %    T= c_mdl.mk_coarse_fine_mapping.f2c_offset (1xN_dims)
0024 %    M= c_mdl.mk_coarse_fine_mapping.f2c_project (N_dimsxN_dims)
0025 %  by default T= [0,0,0] and M=eye(3)
0026 %
0027 % if c_mdl is 2D and f_mdl is 3D, then parameter
0028 %     c_mdl.mk_coarse_fine_mapping.z_depth
0029 %     indicates the +/- z_depth which elements in 2D are
0030 %     considered to be extruded in 3D (default inf)
0031 %
0032 % NOTES:
0033 % if c_mdl and f_mdl do not cover the same area, then
0034 %    sum(c2f') will not be 1. If all coarse elements cover
0035 %    at least partially the fine ones, then this
0036 %    can be corrected by:
0037 %      c2f = c2f./(sum(c2f,2) * ones(1,size(c2f,2)));
0038 % if not all coarse elements cover fine ones, then this
0039 %    approach cannot be used. This will be fixed in a
0040 %    future release
0041 %
0042 % WARNING:
0043 % If c_mdl is not simply connected, the results are wrong!
0044 %
0045 % See also MK_C2F_CIRC_MAPPING
0046 
0047 % (C) 2007-2012 Andy Adler. License: GPL version 2 or version 3
0048 % $Id: mk_approx_c2f.m 5880 2018-12-21 23:19:23Z aadler $
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 % Mapping depends only on nodes and elems - remove the other stuff
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 % find all elems of ff_mdl completely contained in cc_mdl
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     % shrink pts slightly (s_fac) so they're not on the boundary
0092     % by shrinking, we avoid cases where an element is
0093     % only slighly intersecting another. This is beyond the
0094     % resolution of the next step (interpolation) anyway
0095     s_fac= .9; % .9999;
0096     for dim= 1:ef-1 % ef-1 is dimensions in fm
0097        % fm_pts is local_nodes x elems x xyz
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     % if all points are outside (NaN) then c_elems = -1
0107     % if all points are in one elem   then c_elems = elem #
0108     % if all points are in diff elems then c_elems = 0
0109     c_elems= all(diff(tsn,1,2)==0,2) .* tsn(:,1);
0110     c_elems(all(tsn==-1,2))= -1; % all points outside
0111 
0112 
0113 % tsn = vector of length size(fm_pts,1)
0114 % tsn(i) = elem in cm which contains point
0115 % tsn(i) = -1 if point is outside cm (and z_depth, if appropriate)
0116 function tsn= search_fm_pts_in_cm(cm, fm_pts, z_depth);
0117     dc= size(cm.elems,2)-1;  %coarse dim
0118     df= size(fm_pts,2); %fine dim
0119 
0120     tsn= -ones(size(fm_pts,1),1);
0121     not_oor= (tsn==-1); % logical 1
0122 
0123     if dc==2  %corse model is 2D
0124 
0125        if df==3
0126        % look for f_mdl z not out of range
0127           not_oor= not_oor &  any( abs(fm_pts(:,3) ) <= z_depth , 2);
0128        end
0129        dims=1:2;
0130 
0131     elseif dc==3 %coarse model is 3D
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 % find fraction of elem contained in cm
0143 % idx is the index into which the elem is contained
0144 % if idx >= 1, the element is completely with in that coarse elem
0145 % if idx == 0, the element is crosses several coarse elems
0146 % if idx ==-1, the element is outside the coarse model
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    % lower density interpolation in higher dimentions, since
0157    % the added dimensions will give extra interpolation points.
0158    n_interp = 7-df;
0159    interp_mdl.type= 'fwd_model';
0160    interp_mdl.nodes= fm.nodes;
0161    interp_mdl.elems= fm.elems(fidx,:);
0162    interp_mdl.interp_mesh.n_interp = fm.interp_mesh.n_interp;
0163  
0164 
0165    fm_pts = interp_mesh( interp_mdl, n_interp);
0166    l_interp = size(fm_pts,3);
0167 
0168    fm_pts = permute( fm_pts, [3,1,2]);
0169    fm_pts = reshape(fm_pts, [], df-1);
0170 
0171    tsn= search_fm_pts_in_cm(cm, fm_pts, z_depth);
0172 
0173    tsn_idx= ones(l_interp,1)*fidx(:)';
0174    tsn_idx= tsn_idx(:);
0175    % find and isolate outside elements
0176    outside_idx= tsn==-1;
0177    tsn(outside_idx) = [];
0178    tsn_idx(outside_idx) = [];
0179    
0180    in_idx= find(idx<=0);
0181    ridx= 1:nf; ridx(in_idx)=[];
0182    idx(in_idx)=[];
0183 
0184    % first term is contribution from f_elems in one c_elem
0185    % next term is contribution from f_elems in many c_elems, weighted
0186    c_elems = sparse(ridx,idx,1,nf,nc) +  ...
0187              sparse(tsn_idx,tsn,1,nf,nc)/l_interp;
0188       
0189 % Do 3D interpolation of region xyzmin= [x,y,z] to xyzmax
0190 %  with n_interp points in the minimum direction
0191 function xyz = interpxyz( xyzmin, xyzmax, n_interp)
0192     xyzdelta= xyzmax - xyzmin;
0193     xyz_interp = 1 + floor(n_interp * xyzdelta / max(xyzdelta) );
0194     xspace = linspace(xyzmin(1), xyzmax(1), xyz_interp(1) );
0195     yspace = linspace(xyzmin(2), xyzmax(2), xyz_interp(2) );
0196     zspace = linspace(xyzmin(3), xyzmax(3), xyz_interp(3) );
0197     [xx3,yy3,zz3] = ndgrid( xspace, yspace, zspace );
0198     xyz= [xx3(:), yy3(:), zz3(:)];
0199 
0200 % Offset and project f_mdl as required
0201 function f_mdl= offset_and_project( f_mdl, c_mdl)
0202     [fn,fd]= size(f_mdl.nodes);
0203     T= c_mdl.mk_coarse_fine_mapping.f2c_offset;
0204     M= c_mdl.mk_coarse_fine_mapping.f2c_project;
0205     
0206     f_mdl.nodes= (M*( f_mdl.nodes - ones(fn,1)*T )')';
0207 
0208 function [c_mdl f_mdl] = assign_defaults( c_mdl, f_mdl )
0209     [fn,fd]= size(f_mdl.nodes);
0210     try    
0211         c_mdl.mk_coarse_fine_mapping.f2c_offset; % test exist
0212     catch
0213         c_mdl.mk_coarse_fine_mapping.f2c_offset= zeros(1,fd);
0214     end
0215     try    
0216         c_mdl.mk_coarse_fine_mapping.f2c_project;
0217     catch
0218         c_mdl.mk_coarse_fine_mapping.f2c_project= speye(fd);
0219     end
0220     try    
0221         c_mdl.mk_coarse_fine_mapping.z_depth;
0222     catch
0223         c_mdl.mk_coarse_fine_mapping.z_depth= inf;
0224     end
0225     try    
0226         f_mdl.interp_mesh.n_interp;
0227     % lower density interpolation in higher dimentions, since
0228     % the added dimensions will give extra interpolation points.
0229     catch
0230         f_mdl.interp_mesh.n_interp = 7 - size(f_mdl.elems,2);
0231     end
0232      
0233     
0234 function do_unit_test
0235     fmdl = mk_circ_tank(2,[],2); fmdl.nodes = fmdl.nodes*2;
0236     cmdl = mk_circ_tank(2,[],2); cmdl.nodes = cmdl.nodes*2;
0237     c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0238     unit_test_cmp('t1',c2f,eye(16),1e-15)
0239 
0240     fmdl = mk_circ_tank(3,[],2);
0241     fmdl.nodes = fmdl.nodes*3;
0242     c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0243     unit_test_cmp('t2 analytic',c2f,[eye(16);zeros(20,16)],1e-15)
0244     c2f = mk_approx_c2f( fmdl, cmdl);
0245     unit_test_cmp('t2 approx',c2f,[eye(16);zeros(20,16)],1e-15)
0246 
0247     fmdl = mk_circ_tank(2,[],2); fmdl.nodes = fmdl.nodes*2;
0248     cmdl = mk_circ_tank(1,[],2); cmdl.nodes = cmdl.nodes*1;
0249     c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0250     unit_test_cmp('t3 analytic',c2f,[eye(4);zeros(12,4)],1e-15)
0251     c2f = mk_approx_c2f( fmdl, cmdl);
0252     unit_test_cmp('t3 approx',c2f,[eye(4);zeros(12,4)],1e-15)
0253 
0254     fac = 0.8;
0255     cmdl = mk_circ_tank(1,[],2); cmdl.nodes = cmdl.nodes*fac;
0256     c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0257     unit_test_cmp('t4 analytic',c2f,[eye(4)*fac^2;zeros(12,4)],1e-15)
0258     c2f = mk_approx_c2f( fmdl, cmdl);
0259     unit_test_cmp('t4 approx',c2f,[eye(4)*2/3;zeros(12,4)],1e-15)
0260 
0261     fac=1.2;
0262     cmdl = mk_circ_tank(1,[],2); cmdl.nodes = cmdl.nodes*fac;
0263     c2f = mk_approx_c2f( fmdl, cmdl);
0264     unit_test_cmp('t5',c2f,[eye(4);eye(4)/3;kron(eye(4),[1;1])/15],1e-15);
0265 
0266     fmdl = mk_circ_tank(10,[],2);
0267     cmdl = mk_circ_tank(8,[],2);
0268     c2f = mk_approx_c2f( fmdl, cmdl);
0269     unit_test_cmp('t6',sum(c2f'),ones(1,size(c2f,1)),1e-14);
0270    
0271     cmdl.nodes = cmdl.nodes*0.95;
0272 % show_fem(fmdl); hold on ; show_fem(cmdl); hold off
0273     c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0274 
0275 function load
0276 
0277 % Create forward, fine tank model
0278 electrodes_per_plane = 16;
0279 number_of_planes = 2;
0280 tank_radius = 0.2;
0281 tank_height = 0.5;
0282 fine_mdl = ng_mk_cyl_models([tank_height,tank_radius],...
0283     [electrodes_per_plane,0.15,0.35],[0.01]);
0284  
0285 % Create coarse model for inverse problem
0286 coarse_mdl_maxh = 0.07; % maximum element size
0287 coarse_mdl = ng_mk_cyl_models([tank_height,tank_radius,coarse_mdl_maxh],[0],[]);
0288 
0289 disp('Calculating coarse2fine mapping ...');
0290 inv3d.fwd_model.coarse2fine = ...
0291        mk_coarse_fine_mapping( fine_mdl, coarse_mdl);
0292 disp('   ... done');