Home > eidors > models > mk_approx_c2f.m

mk_approx_c2f

PURPOSE ^

MK_APPROX_C2F: create a mapping matrix from coarse to fine FEM

SYNOPSIS ^

function [mapping, outside] = mk_approx_c2f( f_mdl, c_mdl );

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 5599 2017-06-22 02:29:36Z 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.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    % find and isolate outside elements
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    % first term is contribution from f_elems in one c_elem
0184    % next term is contribution from f_elems in many c_elems, weighted
0185    c_elems = sparse(ridx,idx,1,nf,nc) +  ...
0186              sparse(tsn_idx,tsn,1,nf,nc)/l_interp;
0187       
0188 % Do 3D interpolation of region xyzmin= [x,y,z] to xyzmax
0189 %  with n_interp points in the minimum direction
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 % Offset and project f_mdl as required
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; % test exist
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     % lower density interpolation in higher dimentions, since
0227     % the added dimensions will give extra interpolation points.
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 % show_fem(fmdl); hold on ; show_fem(cmdl); hold off
0272     c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0273 
0274 function load
0275 
0276 % Create forward, fine tank model
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 % Create coarse model for inverse problem
0285 coarse_mdl_maxh = 0.07; % maximum element size
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');

Generated on Fri 01-Jun-2018 15:59:55 by m2html © 2005