mk_analytic_c2f

 MK_ANALYTIC_C2F: create a mapping matrix from coarse to fine FEM
 [c2f,out]= mk_analytic_c2f( f_mdl, c_mdl, opt );
  
 Parameters:
    c_mdl is coarse fwd_model
    f_mdl is fine fwd_model
    opt   is option structure

 C2F_ij is the fraction of 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_analytic_c2f.f2c_offset (1xN_dims)
    M= c_mdl.mk_analytic_c2f.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_analytic_c2f.z_depth
     indicates the +/- z_depth which elements in 2D are
     considered to be extruded in 3D (default inf)

 See also MK_C2F_CIRC_MAPPING, MK_APPROX_C2F, MK_COARSE_FINE_MAPPING

0001 function [mapping, outside] = mk_analytic_c2f( f_mdl, c_mdl, opt)
0002 % MK_ANALYTIC_C2F: create a mapping matrix from coarse to fine FEM
0003 % [c2f,out]= mk_analytic_c2f( f_mdl, c_mdl, opt );
0004 %
0005 % Parameters:
0006 %    c_mdl is coarse fwd_model
0007 %    f_mdl is fine fwd_model
0008 %    opt   is option structure
0009 %
0010 % C2F_ij is the fraction of f_mdl element i which is
0011 %   contained in c_mdl element j. This is used to map
0012 %   from data on the reconstruction model (c_mdl) to
0013 %   the forward model f_mdl as
0014 %      elem_data_fine = Mapping*elem_data_coase
0015 %
0016 % OUT_i is the fraction of f_mdl element i which is not
0017 %   contained in any c_mdl element.
0018 %
0019 % OPTIONS:
0020 % if the geometry of the fine and coarse models are not
0021 %  aligned, then they can be translated and mapped using
0022 %    coarse_xyz = (M* (fine_xyz - T)')'
0023 %  where
0024 %    T= c_mdl.mk_analytic_c2f.f2c_offset (1xN_dims)
0025 %    M= c_mdl.mk_analytic_c2f.f2c_project (N_dimsxN_dims)
0026 %  by default T= [0,0,0] and M=eye(3)
0027 %
0028 % if c_mdl is 2D and f_mdl is 3D, then parameter
0029 %     c_mdl.mk_analytic_c2f.z_depth
0030 %     indicates the +/- z_depth which elements in 2D are
0031 %     considered to be extruded in 3D (default inf)
0032 %
0033 % See also MK_C2F_CIRC_MAPPING, MK_APPROX_C2F, MK_COARSE_FINE_MAPPING
0034 
0035 % (C) 2007-2015 Andy Adler and Bartlomiej Grychtol.
0036 % License: GPL version 2 or version 3
0037 % $Id: mk_analytic_c2f.m 5884 2018-12-22 14:16:36Z aadler $
0038 
0039 if ischar(f_mdl) && strcmp(f_mdl, 'UNIT_TEST'); do_unit_test; return; end
0040 
0041 if ischar(f_mdl) && strcmp(f_mdl, 'LOAD'); preload; return; end
0042 
0043 if nargin== 2, opt = struct; end
0044 
0045 opt = assign_defaults( c_mdl, f_mdl, opt );
0046 
0047 copt.cache_obj = cache_obj(c_mdl, f_mdl, opt);
0048 copt.fstr = 'mk_analytic_c2f';
0049 
0050 [mapping, outside] = eidors_cache(@mapping_calc,{f_mdl,c_mdl,opt},copt);
0051 
0052 
0053 
0054 function [mapping, outside] = mapping_calc(f_mdl, c_mdl, opt)
0055     
0056    f_mdl= offset_and_project( f_mdl, opt);
0057  
0058    f_dim = elem_dim(f_mdl);
0059    c_dim = elem_dim(c_mdl);
0060    d_num = f_dim * 10 + c_dim;
0061    
0062    switch d_num
0063       case 33
0064          mapping = mk_tet_c2f(f_mdl,c_mdl, opt);
0065       case 22
0066          mapping = mk_tri_c2f(f_mdl,c_mdl, opt);
0067       case 32
0068          mapping = mk_tri2tet_c2f(f_mdl, c_mdl, opt);
0069       case 23
0070          mapping = mk_tri2tet_c2f(c_mdl, f_mdl, opt);
0071          mapping = bsxfun(@times,   mapping , get_elem_volume(c_mdl));
0072          if isinf(opt.z_depth)
0073             height = max(c_mdl.nodes(:,3)) - min(c_mdl.nodes(:,3));
0074          else
0075             height = opt.z_depth;
0076          end
0077          mapping = mapping * height;
0078          mapping = bsxfun(@rdivide, mapping', get_elem_volume(f_mdl)); 
0079       otherwise
0080          error('EIDORS:WrongDim',['Unsupported combination of dimensions: ' ...
0081             'c_mdl=%d, f_mdl=%d'],c_dim,f_dim);
0082    end
0083 
0084    
0085    if isfield(c_mdl,'coarse2fine')
0086       mapping = mapping*c_mdl.coarse2fine;
0087    end
0088    
0089    if nargout>1;
0090       outside = 1 - sum(mapping,2);
0091    end
0092 
0093 
0094 % Mapping depends only on nodes and elems - remove the other stuff
0095 function c_obj = cache_obj(c_mdl, f_mdl, opt)
0096    c_obj = {c_mdl.nodes, c_mdl.elems,  ...
0097             f_mdl.nodes, f_mdl.elems, opt};
0098 
0099 % Offset and project f_mdl as required
0100 function f_mdl= offset_and_project( f_mdl, opt)
0101     [fn,fd]= size(f_mdl.nodes);
0102     T= opt.f2c_offset;
0103     M= opt.f2c_project;
0104     
0105     f_mdl.nodes= (M*( f_mdl.nodes - ones(fn,1)*T )')';
0106 
0107 function opt = assign_defaults( c_mdl, f_mdl, org )
0108    opt = struct;
0109    try opt = c_mdl.mk_analytic_c2f; end
0110    fn = fieldnames(org);
0111    for f = 1:numel(fn)
0112       opt.(fn{f}) = org.(fn{f}); % explicit options override those on the model
0113    end
0114    [fn,fd]= size(f_mdl.nodes);
0115    try    opt.f2c_offset; % test exist
0116    catch  opt.f2c_offset= zeros(1,fd);
0117    end
0118    try    opt.f2c_project;
0119    catch  opt.f2c_project= speye(fd);
0120    end
0121    try    opt.z_depth;
0122    catch  opt.z_depth= inf;
0123    end
0124 
0125      
0126     
0127 function do_unit_test
0128    ll = eidors_msg('log_level', 1);
0129    do_original_2d_tests
0130    test_2d3d_handling
0131    mk_tri_c2f UNIT_TEST
0132    mk_tet_c2f UNIT_TEST
0133    mk_tri2tet_c2f UNIT_TEST
0134    eidors_msg('log_level',ll);
0135    
0136 function test_2d3d_handling
0137    tet = mk_grid_model([],[0 1],[0 1],[0 1]);
0138    tet = rmfield(tet,'coarse2fine');
0139    tri = mk_grid_model([],[0 1],[0 1]);
0140    tri = rmfield(tri,'coarse2fine');
0141    
0142    c2f_a = mk_analytic_c2f(tri,tet);
0143    c2f_a = bsxfun(@times, c2f_a, get_elem_volume(tri));
0144    c2f_b = mk_analytic_c2f(tet,tri);
0145    c2f_b = bsxfun(@times, c2f_b, get_elem_volume(tet));
0146    
0147    unit_test_cmp('2d3d:01',sum(c2f_a(:)), 1,eps);
0148    unit_test_cmp('2d3d:02',sum(c2f_b(:)), 1,eps);
0149    unit_test_cmp('2d3d:03',c2f_a,c2f_b',eps)
0150    
0151 
0152    
0153 function do_original_2d_tests
0154    fmdl = mk_circ_tank(2,[],2); fmdl.nodes = fmdl.nodes*2;
0155    cmdl = mk_circ_tank(2,[],2); cmdl.nodes = cmdl.nodes*2;
0156    c2f = mk_analytic_c2f( fmdl, cmdl);
0157    unit_test_cmp('t1',c2f,eye(16),eps)
0158 
0159    fmdl = mk_circ_tank(3,[],2);
0160    fmdl.nodes = fmdl.nodes*3;
0161    c2f = mk_analytic_c2f( fmdl, cmdl);
0162    unit_test_cmp('t2',c2f,[eye(16);zeros(20,16)],eps)
0163 
0164    fmdl = mk_circ_tank(2,[],2); fmdl.nodes = fmdl.nodes*2;
0165    cmdl = mk_circ_tank(1,[],2); cmdl.nodes = cmdl.nodes*1;
0166    c2f = mk_analytic_c2f( fmdl, cmdl);
0167    unit_test_cmp('t3',c2f,[eye(4);zeros(12,4)],eps)
0168 
0169    cmdl = mk_circ_tank(1,[],2); cmdl.nodes = cmdl.nodes*0.8;
0170    c2f = mk_analytic_c2f( fmdl, cmdl);
0171    unit_test_cmp('t4',c2f,[eye(4)*.64;zeros(12,4)], eps)
0172 
0173 
0174    fmdl = mk_circ_tank(10,[],2);
0175    cmdl = mk_circ_tank(8,[],2);
0176    c2f = mk_analytic_c2f( fmdl, cmdl);
0177    unit_test_cmp('t6',sum(c2f(1:324,:)'),ones(1,324),1e-14);
0178 
0179    cmdl.nodes = cmdl.nodes*0.95;
0180    % show_fem(fmdl); hold on ; show_fem(cmdl); hold off
0181    c2f = mk_analytic_c2f( fmdl, cmdl);
0182 
0183 % preload into cache before tests
0184 function preload
0185 
0186    % Create forward, fine tank model
0187    electrodes_per_plane = 16;
0188    number_of_planes = 2;
0189    tank_radius = 0.2;
0190    tank_height = 0.5;
0191    fine_mdl = ng_mk_cyl_models([tank_height,tank_radius],...
0192        [electrodes_per_plane,0.15,0.35],[0.01]);
0193     
0194    % Create coarse model for inverse problem
0195    coarse_mdl_maxh = 0.07; % maximum element size
0196    coarse_mdl = ng_mk_cyl_models([tank_height,tank_radius,coarse_mdl_maxh],[0],[]);
0197 
0198    disp('Calculating coarse2fine mapping ...');
0199    inv3d.fwd_model.coarse2fine = ...
0200           mk_analytic_c2f( fine_mdl, coarse_mdl);
0201    disp('   ... done');