mk_c2f_circ_mapping

 MK_C2F_CIRC_MAPPING: create a mapping matrix from circles/spheres to FEM
 mapping= mk_c2f_circ_mapping( mdl, xyzr );

 Mapping approximates elem_data_fine from elem_data_coase
   elem_data_model = Mapping * elem_data_circles

 mdl is coarse fwd_model
 xyzr is the 3xN matrix (2D) or 4xN matrix (3D) of
      circle centres and radii

 this function approximates using points interpolated into elements
   use mdl.interp_mesh.n_points to control interpolation density  

 if a 3xN matrix is specified for a 3D model, then cylindrical
  shapes (circle extruded in z) are selected

0001 function [mapping failed] = mk_c2f_circ_mapping( mdl, xyzr );
0002 % MK_C2F_CIRC_MAPPING: create a mapping matrix from circles/spheres to FEM
0003 % mapping= mk_c2f_circ_mapping( mdl, xyzr );
0004 %
0005 % Mapping approximates elem_data_fine from elem_data_coase
0006 %   elem_data_model = Mapping * elem_data_circles
0007 %
0008 % mdl is coarse fwd_model
0009 % xyzr is the 3xN matrix (2D) or 4xN matrix (3D) of
0010 %      circle centres and radii
0011 %
0012 % this function approximates using points interpolated into elements
0013 %   use mdl.interp_mesh.n_points to control interpolation density
0014 %
0015 % if a 3xN matrix is specified for a 3D model, then cylindrical
0016 %  shapes (circle extruded in z) are selected
0017 %
0018 
0019 % (C) 2009 Andy Adler. License: GPL version 2 or version 3
0020 % $Id: mk_c2f_circ_mapping.html 2819 2011-09-07 16:43:11Z aadler $
0021 
0022 if isstr(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0023 
0024 failed = false;% not all subfunctions set this
0025 
0026 c_obj = cache_obj(mdl, xyzr);
0027 
0028 mapping = eidors_obj('get-cache', c_obj, 'circle_mapping');
0029 failed = eidors_obj('get-cache', c_obj, 'failed_circle_mapping');
0030 if ~isempty(mapping)
0031     eidors_msg('mk_c2f_circ_mapping: using cached value', 3);
0032 else
0033 
0034     mdl = fix_model(mdl);
0035     switch size(xyzr,1)
0036       case 3; mapping = contained_elems_2d( mdl, xyzr );
0037       case 4; [mapping failed] = contained_elems_3d( mdl, xyzr );
0038       otherwise; error('size of xyzr incorrect');
0039     end
0040 
0041     eidors_obj('set-cache', c_obj, 'circle_mapping', mapping);
0042     eidors_obj('set-cache', c_obj, 'failed_circle_mapping', failed);
0043     eidors_msg('mk_coarse_fine_mapping: setting cached value', 3);
0044 end
0045 
0046 % Mapping depends only on nodes and elems - remove the other stuff
0047 function c_obj = cache_obj(mdl, xyzr)
0048    c_obj = {mdl.nodes, mdl.elems, xyzr};
0049 
0050 % Redirector during test code dev
0051 function mapping = contained_elems_2d( mdl, xyr );
0052 %mapping = contained_elems_2d_new( mdl, xyr );
0053  mapping = contained_elems_2d_old( mdl, xyr );
0054 
0055 function mapping = contained_elems_2d_new( mdl, xyr );
0056    % We fill sparse by columns, (ie adding in CCS storage)
0057    Nc = size(xyr,      2); % Num circs
0058    too_far = elems_too_far( mdl, xyr );
0059 
0060    mapping = sparse( num_elems(mdl) , Nc );
0061    for i=1:Nc
0062      mapping(:,i) = circ_in_elem_2d(mdl, find( ~too_far(:,i)), ...
0063                 xyr(1,i), xyr(2,i), xyr(3,i));
0064    end
0065 
0066 % look only at elements 'look'
0067 function mapping = circ_in_elem_2d( mdl, look, xc, yc, rc);
0068    Nt = elem_dim(mdl) + 1; % nodes per simplex
0069    pirc2 = pi*rc^2;
0070 % start assuming no content
0071    mapping = sparse(num_elems(mdl),1);
0072 % For each element, find nodes inside
0073    els = mdl.elems(look,:);
0074    ndx = reshape(mdl.nodes(els,1) - xc, size(els));
0075    ndy = reshape(mdl.nodes(els,2) - yc, size(els));
0076    n_in = (ndx.^2 + ndy.^2) < rc^2; % node inside
0077 % triangles with 3 nodes inside are all in, stop looking at them
0078    all_n_in = sum(n_in,2) == Nt;
0079    mapping(look(all_n_in)) = 1;
0080    look(all_n_in)= []; n_in(all_n_in,:)= [];
0081 % find distance inside each face
0082    f_in = zeros( length(look), Nt); 
0083    k=1;  for i= look(:)';
0084       faces = mdl.elem2face(i,:);
0085       out   =~mdl.inner_normal(i,:);
0086       f_norm= mdl.normals( faces, :);
0087       f_norm(out,:) = -f_norm(out,:);
0088 
0089       f_ctr = mdl.face_centre( faces,:);
0090       v_ctr = repmat([xc,yc],Nt,1) - f_ctr;
0091       v_ctr = sum(v_ctr .* f_norm,2)/rc; % >1 in, <-1 out,
0092       f_in(k,:) = v_ctr';
0093    k=k+1;end
0094 
0095 % triangles with any sides outside are out
0096    any_s_out= any(f_in<-1,2);
0097    look(any_s_out)= [];
0098    n_in(any_s_out,:) = [];
0099    f_in(any_s_out,:) = [];
0100 
0101 % triangles with 3 sides inside are all in.
0102    all_s_in = sum(f_in>1,2) == Nt;
0103    mapping(look(all_s_in)) = pirc2 / ...
0104             mdl.elem_volume(look(all_s_in));
0105    look(all_s_in)= [];
0106    n_in(all_s_in) = [];
0107    f_in(all_s_in) = [];
0108 
0109 % Now, all triangles should be partially in
0110 % Calculate area chopped out
0111    fin1 = f_in<1;
0112    a_out = zeros(size(fin1));
0113    a_out(fin1) = acos(f_in(fin1));
0114    a_out(fin1) = (a_out(fin1) - cos(a_out(fin1)).*sin(a_out(fin1)))/pi;
0115    
0116    % start with the default. This is accurate if there are
0117    % no contained nodes, otherwise we need to add back
0118    % those fractions
0119    mapping(look) = pirc2 / mdl.elem_volume(look);
0120 
0121    %TODO: rewrite loop to avoid case 0.
0122    k=1; for i= look(:)';
0123       vol = pi*rc^2 / mdl.elem_volume(i);
0124       switch sum(n_in(k,:))
0125          case 0; % already do this
0126    
0127          case 1; 
0128             nd = mdl.elems(k, n_in(k,:));
0129    keyboard
0130             vol = vol + pi_slice(p1,p2,[xc,yc],mdl.nodes(nd,:),rc);
0131 
0132          case 2; 
0133             nd = mdl.elems(k, n_in(k,:));
0134             vol = vol ...
0135              + pi_slice(p1,p2,[xc,yc],mdl.nodes(nd(1),:),rc) ...
0136              + pi_slice(p1,p2,[xc,yc],mdl.nodes(nd(2),:),rc);
0137 
0138          otherwise; error('cant get here'); 
0139       end
0140    k=k+1; end
0141    
0142 
0143 % Calculate the area of a slice
0144 function a = pi_slice(p1,p2,c,p,r)
0145   a_p12 = 0.5*abs(det([1,p1;1,p2;1,p]));
0146 
0147   a_c12 = 0.5*abs(det([1,p1;1,p2;1,c]));
0148   np1c  = p1-c; np1c = np1c / norm(np1c);
0149   np2c  = p2-c; np2c = np2c / norm(np2c);
0150   ang   = acos( dot(np1c,np2c) );
0151   area  = ang*r^2/2 - a_c12 + a_p12;
0152      
0153 
0154 function mapping = contained_elems_2d_old( mdl, xyr );
0155    Ne = size(mdl.elems,1); % Num elems
0156    Nc = size(xyr,      2); % Num circs
0157    % We fill sparse by columns, due to CCS storage, this is fairly efficient
0158    mapping = sparse( Ne, Nc );
0159 
0160    % Interpolate
0161    n_interp =  7-size(mdl.nodes,2);
0162    m_pts = interp_mesh( mdl, n_interp); 
0163    for i=1:Nc
0164      xc = m_pts(:,1,:) - xyr(1,i);
0165      yc = m_pts(:,2,:) - xyr(2,i);
0166      inr= xc.^2 + yc.^2 < xyr(3,i)^2;
0167      frac= mean(inr,3);
0168      mapping(:,i) = frac;
0169    end
0170 
0171    % 1. Get furthest node in each element
0172    % 2. Get the max edge length of each elem
0173    % 3. Find elems that are too far
0174 function too_far = elems_too_far( mdl, xyr );
0175    Ne = num_elems(mdl);
0176    Nc = size(xyr, 2); % Num circs
0177    Nt = elem_dim(mdl) + 1; % Elements per simplex
0178    nodes = repmat(mdl.nodes,[1,1,Nc]);
0179    targets = repmat(xyr(1:mdl_dim(mdl),:),[1,1,num_nodes(mdl)]);
0180    targets = shiftdim(targets,2);
0181    dist = nodes - targets;
0182    dist = sqrt(sum(dist.^2,2));
0183    node_target_dist = squeeze(dist);
0184    furthest_elem_node_dist = node_target_dist(mdl.elems,:);
0185    furthest_elem_node_dist = reshape(furthest_elem_node_dist,Ne,Nt,Nc);
0186    [furthest_elem_node_dist, furthest_elem_node]= max(furthest_elem_node_dist,[],2);
0187    furthest_elem_node_dist = squeeze(furthest_elem_node_dist);
0188    furthest_elem_node = squeeze(furthest_elem_node);
0189    
0190    max_edge_len =  repmat(mdl.max_edge_len,1,Nc);
0191    radius = ones(Ne,1)*xyr(Nt,:);
0192    too_far = (furthest_elem_node_dist - max_edge_len) > radius;
0193 
0194 function [mapping failed] = contained_elems_3d( mdl, xyr );
0195    Ne = size(mdl.elems,1); % Num elems
0196    Nc = size(xyr,      2); % Num circs
0197    failed(1:Nc) = false;
0198    % We fill sparse by columns, due to CCS storage, this is fairly efficient
0199    mapping = sparse( Ne, Nc );
0200     if 0
0201    % INterpolate
0202    n_interp = 4; % 7-df
0203    m_pts = interp_mesh( mdl, n_interp); 
0204    for i=1:Nc
0205      mapping(:,i) = contained_elem_pts(m_pts, xyr(:,i));
0206    end
0207     else
0208 
0209    % 4. Make a tmp model with only the remaining elems
0210    % 5. Interpolate
0211    % 6. Merge
0212 
0213    too_far = elems_too_far( mdl, xyr );
0214    
0215    tmp = eidors_obj('fwd_model','tmp','nodes',mdl.nodes,'elems',mdl.elems);
0216    %mapping = sparse( Ne, Nc );
0217    n_interp_min = 6;
0218    n_interp_max = 10;
0219    for i=1:Nc
0220        good = ~too_far(:,i);
0221        if ~any(good), continue, end %point outside the mesh
0222        tmp.elems = mdl.elems(good,:);
0223        n_interp = n_interp_min-1;
0224        log_level = eidors_msg('log_level',1);
0225        while(sum(mapping(good,i))==0 && n_interp < n_interp_max-1)
0226            n_interp = n_interp+1;
0227            m_pts = interp_mesh( tmp, n_interp);
0228            mapping(good,i) = contained_elem_pts(m_pts, xyr(:,i));
0229        end
0230        eidors_msg('log_level', log_level);
0231        if (sum(mapping(good,i)) == 0)
0232            failed(i) = true;
0233            eidors_msg(['mk_c2f_circ_mapping: Interpolation failed for point ' num2str(i)]);
0234        end
0235    end
0236     end
0237    
0238    
0239 function frac= contained_elem_pts(m_pts, xyr);
0240 % This is more clear
0241 %    xc = m_pts(:,1,:) - xyr(1);
0242 %    yc = m_pts(:,2,:) - xyr(2);
0243 %    zc = m_pts(:,3,:) - xyr(3);
0244 %    inr= xc.^2 + yc.^2 + zc.^2 < xyr(4)^2;
0245 
0246 % But this is how to stop matlab from wasting memory
0247      inr = (m_pts(:,1,:) - xyr(1)).^2 + ...% xc =
0248            (m_pts(:,2,:) - xyr(2)).^2 + ...% yc =
0249            (m_pts(:,3,:) - xyr(3)).^2;     % zc =
0250      inpts = inr < xyr(4)^2;
0251 
0252 %    frac= mean( inpts ,3);
0253 %    FIXME: Octave doesn't like to mean on logical
0254      frac= mean( int8( inpts ) ,3);
0255      if sum(inpts(:))==0
0256          % TODO: This message is outdated
0257          eidors_msg(['mk_c2f_circ_mapping: Interpolation failed: increase ', ...
0258                          'fwd_model.interp_mesh.n_interp']);
0259      end
0260 
0261 function do_unit_test
0262    %2D example
0263    imdl = mk_common_model('a2c2',16); fmdl=imdl.fwd_model;
0264    xyc = [0,0.27,0.18;0,-0.1,0.03;0,0.1,0.2;0.1,0.37,0.1]';
0265    th=linspace(0,2*pi,20)';
0266    xx=[0*th+1]*xyc(1,:)+sin(th)*xyc(3,:);
0267    yy=[0*th+1]*xyc(2,:)+cos(th)*xyc(3,:);
0268    show_fem(fmdl,[0,0,1]); set(line(xx,yy),'LineWidth',2);
0269 
0270    c2f= mk_c2f_circ_mapping( fmdl, [0;0;0.1] );
0271    t1= all( abs(c2f(1:4)-0.2857)<.001 ) & all( c2f(5:end)==0 );
0272    unit_test_cmp('2D ex 1:',t1,1);
0273 
0274    c2f= mk_c2f_circ_mapping( fmdl, [.0;.05;0.03]); 
0275    t2= abs( c2f(1) - 0.1429) < .001 & all( c2f(2:end)==0 );
0276    unit_test_cmp('2D ex 2:',t2,1);
0277       
0278    %3D example - cylinder
0279    imdl = mk_common_model('a3cr',16); fmdl=imdl.fwd_model;
0280    c2f= mk_c2f_circ_mapping( fmdl, [0;0;0.1]); 
0281    t3= all( abs(c2f(1:4)-0.1714)<.001 ) & all( c2f(5:64)==0 );
0282    unit_test_cmp('3D ex 1:',t2,1);
0283 
0284    %3D example - cylinder
0285    imdl = mk_common_model('a3cr',16); fmdl=imdl.fwd_model;
0286    c2f= mk_c2f_circ_mapping( fmdl, [0;0;0;0.1]); 
0287    t4= all( abs(c2f(193:196)-0.0595)<.001 ) & all( c2f(1:64)==0 );
0288    unit_test_cmp('3D ex 2:',t4,1);
0289    
0290    %3D example - cylinder - 2 pts
0291    imdl = mk_common_model('a3cr',16); fmdl=imdl.fwd_model;
0292    c2f= mk_c2f_circ_mapping( fmdl, [0 0;0 0;0 0;0.1 0.2]); 
0293    t4= all( abs(c2f(193:196,1)-0.0595)<.001 ) & all( c2f(1:64)==0 );
0294    unit_test_cmp('3D ex 3:',t4,1);