calc_model_reduction

 calc_model_reduction: calculate the fields for a reduced model
    which should speed up forward calculations
 
 mr = calc_model_reducton(fmdl)
  where
    mr.main_region = vector, and 
    mr.regions = struct

 Model Reduction: use precomputed fields to reduce the size of
    the forward solution. Nodes which are 1) not used in the output
    (i.e. not electrodes) 2) all connected to the same conductivity via
    the c2f mapping are applicable.

0001 function mr = calc_model_reduction(fmdl)
0002 % calc_model_reduction: calculate the fields for a reduced model
0003 %    which should speed up forward calculations
0004 %
0005 % mr = calc_model_reducton(fmdl)
0006 %  where
0007 %    mr.main_region = vector, and
0008 %    mr.regions = struct
0009 %
0010 % Model Reduction: use precomputed fields to reduce the size of
0011 %    the forward solution. Nodes which are 1) not used in the output
0012 %    (i.e. not electrodes) 2) all connected to the same conductivity via
0013 %    the c2f mapping are applicable.
0014 
0015 % see: Model Reduction for FEM Forward Solutions, Adler & Lionheart, EIT2016
0016 
0017 % $Id: calc_model_reduction.m 5789 2018-05-21 23:39:41Z aadler $
0018 
0019 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0020 
0021 switch fmdl(1).type
0022   case 'image';      fmdl = fmdl.fwd_model;
0023   case 'inv_model';  fmdl = fmdl.fwd_model;
0024   case 'fwd_model';  fmdl = fmdl;
0025   otherwise;
0026       error('cant process model of type %s', fmdl.type );
0027 end
0028 
0029 if ~isfield(fmdl,'coarse2fine');
0030    FC = system_mat_fields(fmdl);
0031    mr.main_region = logical(ones(1,size(FC,2))'); 
0032    mr.regions     = struct([]);
0033    return
0034 end % not needed without
0035 
0036 
0037 copt.cache_obj.elems = fmdl.elems;
0038 copt.cache_obj.nodes = fmdl.nodes;
0039 copt.cache_obj.coarse2fine = fmdl.coarse2fine;
0040 copt.fstr = 'calc_model_reducton';
0041 copt.log_level = 4;
0042 mr= eidors_cache(@calc_model_reduction_fn,{fmdl},copt );
0043 
0044 function mr = calc_model_reduction_fn(fwd_model);
0045    ne = num_elems(fwd_model);
0046    el = fwd_model.elems; ei=(1:ne)';
0047    e2n = sparse(el,ei*ones(1,elem_dim(fwd_model)+1),1);
0048 
0049    c2n = e2n*fwd_model.coarse2fine;
0050 
0051 % get the number of the c2n mapping. But only work on those with one
0052    ac2n = c2n>0;
0053    ac2n(sum(ac2n,2)>1,:) = 0;
0054    region = ac2n*(1:size(ac2n,2))';
0055    region0 = region==0;
0056    FC = system_mat_fields(fwd_model); S= FC'*FC;
0057    region0(end+1:size(FC,2)) = logical(1);
0058    k=0;
0059    for i=1:max(region);
0060       ff=find(region==i);
0061       if length(ff)>1;
0062          imgn.node_data(ff) = 1+0.5*(k)/5; k=k+1;
0063         %invEi = S(region0,ff)*inv(S(ff,ff))*S(ff,region0);
0064          invEi =(S(region0,ff)/S(ff,ff))*S(ff,region0);
0065          regions(k) = struct('nodes',ff,'field',i,'invE',invEi);
0066       end;
0067    end
0068 
0069    mr.main_region = region0;
0070    mr.regions = regions;
0071 
0072 function do_unit_test
0073    for i=4:4; switch i;
0074       case 1; str = 'a2c0';  tmdl = mk_common_model(str,16);
0075       case 2; str = 'b2C2';  tmdl = mk_common_model(str,16);
0076       case 3; str = 'a3cr';  tmdl = mk_common_model(str,16);
0077       case 4; str = 'ngcyl'; tmdl = ng_mk_cyl_models([1,1,.1],[16,0.5],0.1);
0078    end
0079    img = mk_image(tmdl,1);
0080    img.fwd_model.electrode = img.fwd_model.electrode([15,16,1:3]);
0081    stim = mk_stim_patterns(5,1,[0,1],[0,1],{},1);
0082    img.fwd_model.stimulation = stim;
0083 
0084    pos = interp_mesh(img.fwd_model);
0085    bl = find(pos(:,1)<=0.2 & pos(:,2)<=0.5);
0086    br = find(pos(:,1)>=0.2 & pos(:,2)<=0.5);
0087    img.elem_data(bl) = 0.9;
0088    img.elem_data(br) = 1.1;
0089    pgnode = zeros(1,mdl_dim(img)); pgnode(2) = 0.8;
0090    [~,img.fwd_model.gnd_node] = min(sum( ...
0091                bsxfun(@minus,img.fwd_model.nodes,pgnode).^2,2));
0092    c2f = 1+(1:num_elems(img));
0093    c2f(bl) = 1;
0094    c2f(br) = 2+num_elems(img);
0095    [~,idx2,idx] = unique(c2f); nc = length(idx2);
0096    c2f = sparse(1:num_elems(img),idx,1);
0097 
0098    img.elem_data = c2f*linspace(.7,1.3,length(idx2))';
0099 
0100    vs1 = fwd_solve( img);
0101    
0102 
0103    imgmr = img; imgmr.fwd_model.model_reduction = @calc_model_reduction;
0104    tic;
0105    vs2 = fwd_solve( imgmr);
0106    fprintf('t= %5.3fs:',toc);
0107    unit_test_cmp([str,': 1-2'],vs1.meas,vs2.meas,1e-10);
0108 
0109    img.fwd_model.coarse2fine = c2f;
0110    imgm3 = img;
0111    imgm3.fwd_model.model_reduction = calc_model_reduction( imgm3.fwd_model);
0112    tic;
0113    vs3 = fwd_solve( imgm3);
0114    fprintf('t= %5.3fs:',toc);
0115    unit_test_cmp([str,': 1-3'],vs1.meas,vs3.meas,1e-10);
0116 
0117    imgm4 = img;
0118    imgm4.fwd_model.model_reduction = @calc_model_reduction;
0119    tic;
0120    vs4 = fwd_solve( imgm4);
0121    fprintf('t= %5.3fs:',toc);
0122    unit_test_cmp([str,': 1-4'],vs1.meas,vs3.meas,1e-10);
0123    end