crop_model

 CROP_MODEL: Crop away parts of a fem model

 USAGE #1: crop display to show model internals
   crop_model( axis_handle, fcn_handle );

   fcn_handle ==1 where model is *kept*
 
   if axis_handle==[], then use the current axis
   examples:
     crop_model([],  inline('z==3','x','y','z'))
     crop_model(gca, inline('x+y>0','x','y','z'))

 USAGE #2: crop fwd_model to create new fwd_model
   fmdl_new= crop_model( fwd_model, fcn_handle );
 
   example:
   fmdl2= crop_model(fmdl1, inline('x+y>0','x','y','z'))

  with two parameters output
 [fmdl,c2f_idx]= crop_model( axis_handle, fcn_handle );
     c2f_idx maps each elemen in fmdl_new to fwd_model

0001 function [fmdl,c2f_idx]= crop_model( axis_handle, fcn_handle );
0002 % CROP_MODEL: Crop away parts of a fem model
0003 %
0004 % USAGE #1: crop display to show model internals
0005 %   crop_model( axis_handle, fcn_handle );
0006 %
0007 %   fcn_handle ==1 where model is *kept*
0008 %
0009 %   if axis_handle==[], then use the current axis
0010 %   examples:
0011 %     crop_model([],  inline('z==3','x','y','z'))
0012 %     crop_model(gca, inline('x+y>0','x','y','z'))
0013 %
0014 % USAGE #2: crop fwd_model to create new fwd_model
0015 %   fmdl_new= crop_model( fwd_model, fcn_handle );
0016 %
0017 %   example:
0018 %   fmdl2= crop_model(fmdl1, inline('x+y>0','x','y','z'))
0019 %
0020 %  with two parameters output
0021 % [fmdl,c2f_idx]= crop_model( axis_handle, fcn_handle );
0022 %     c2f_idx maps each elemen in fmdl_new to fwd_model
0023 
0024 % (C) 2006-2008 Andy Adler. License: GPL version 2 or version 3
0025 % $Id: crop_model.m 5112 2015-06-14 13:00:41Z aadler $
0026 
0027 if ischar(axis_handle) && strcmp(axis_handle,'UNIT_TEST'); do_unit_test; return; end
0028 
0029 % TODO (update 2 apr 2012):
0030 %  - make crop_model work for 2D fems
0031 
0032 usage_graphics= 1;
0033 try if axis_handle.type == 'fwd_model'
0034    usage_graphics= 0;
0035 end; end
0036 
0037 if usage_graphics
0038    if isempty(axis_handle)
0039       axis_handle= gca;
0040    end
0041    crop_graphics_model(axis_handle, fcn_handle);
0042 else
0043    [fmdl,c2f_idx]= crop_fwd_model(axis_handle, fcn_handle);
0044 end
0045 
0046 % CROP GRAPHICS
0047 function crop_graphics_model(axis_handle, fcn_handle);
0048    kk= get(axis_handle,'Children');
0049    % only the FEM frame
0050    %k=kk( find( kk== min(kk) ));
0051 
0052    for k= kk(:)'
0053       try
0054          x= get(k,'XData');
0055          y= get(k,'YData');
0056          z= get(k,'ZData');
0057          c= get(k,'CData');
0058          idx= ~all( feval(fcn_handle,x,y,z) );
0059          if any(size(c)>[1,1])
0060             set(k,'Xdata', x(:,idx), ...
0061                   'Ydata', y(:,idx), ...
0062                   'Zdata', z(:,idx), ...
0063                   'Cdata', c(:,idx));
0064          else
0065             set(k,'Xdata', x(:,idx), ...
0066                   'Ydata', y(:,idx), ...
0067                   'Zdata', z(:,idx));
0068          end
0069       end
0070    end
0071 
0072 % CROP fwd_model
0073 function [fmdl1,c2f_idx]= crop_fwd_model(fmdl0, fcn_handle);
0074    fmdl1= fmdl0;
0075 
0076 % Find nodes to remove
0077    nodes= fmdl0.nodes;
0078    [n,d]= size(nodes);
0079    n2xyz= eye(d,3); 
0080    xyz= nodes*n2xyz;
0081 % Keep these nodes
0082    idx0=  all( feval(fcn_handle,xyz(:,1), ...
0083                                 xyz(:,2), ...
0084                                 xyz(:,3)),2);
0085 % remove these nodes
0086    fmdl1.nodes(idx0,:) = [];
0087 
0088 % renumber nodes, set unused ones to 0
0089    idx1= zeros(n,1);
0090    idx1(~idx0)= 1:sum(~idx0);
0091 
0092    fmdl1.elems(:) = idx1(fmdl0.elems);
0093    remove= any( fmdl1.elems == 0, 2);
0094    fmdl1.elems(remove,:)= [];
0095 % c2f_idx maps each new elem to its original position
0096    c2f_idx= find(~remove);
0097 
0098    fmdl1.boundary(:) = idx1(fmdl0.boundary);
0099    remove= any( fmdl1.boundary == 0, 2);
0100    fmdl1.boundary(remove,:)= [];
0101 
0102 % renumber nodes, set unused ones to 0
0103 if isfield(fmdl1,'electrode');
0104    n_elecs = length(fmdl1.electrode);
0105    rm_elec_list = zeros(n_elecs,1);
0106    for i=1:n_elecs;
0107       el_nodes= fmdl0.electrode(i).nodes;
0108       el_nodes(:)= idx1(el_nodes);
0109       if any(el_nodes==0);
0110          rm_elec_list(i) = 1;
0111       end
0112       fmdl1.electrode(i).nodes= el_nodes;
0113    end
0114    if any(rm_elec_list)
0115       str = sprintf('%d,', find(rm_elec_list));
0116       eidors_msg('crop_model: removing electrodes (%s)',str(1:end-1),1);
0117       fmdl1.electrode = fmdl1.electrode( find(~rm_elec_list));
0118    end
0119 end
0120 
0121 
0122 function do_unit_test
0123    subplot(331)
0124    imdl = mk_common_model('n3r2',[16,2]); fmdl= imdl.fwd_model;
0125    show_fem(fmdl);
0126    subplot(332)
0127    show_fem(fmdl); hh= gca; 
0128    subplot(333)
0129    show_fem(fmdl);
0130    crop_model([],inline('z<2','x','y','z'));
0131    crop_model(hh,inline('x<0','x','y','z'));
0132 
0133    subplot(334)
0134    fmdlc = fmdl;
0135    fmdlc = crop_model(fmdlc,inline('x<0','x','y','z'));
0136    show_fem(fmdlc);
0137 
0138    subplot(337)
0139    imdl = mk_common_model('d2c2'); fmdl= imdl.fwd_model;
0140    show_fem(fmdl);
0141    subplot(338)
0142    show_fem(fmdl); hh= gca; 
0143    title('expected fail');
0144    subplot(339)
0145    show_fem(fmdl);
0146    crop_model([],inline('y<0','x','y','z'));
0147    crop_model(hh,inline('x<0','x','y','z'));
0148    title('expected fail');
0149    
0150