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 and Bartlomiej Grychtol.
0020 % License: GPL version 2 or version 3
0021 % $Id: mk_c2f_circ_mapping.m 3934 2013-04-20 15:46:00Z aadler $
0022 
0023 if isstr(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0024 
0025 failed = false;% not all subfunctions set this
0026 
0027 c_obj = cache_obj(mdl, xyzr);
0028 
0029 mapping = eidors_obj('get-cache', c_obj, 'circle_mapping');
0030 failed = eidors_obj('get-cache', c_obj, 'failed_circle_mapping');
0031 if ~isempty(mapping)
0032     eidors_msg('mk_c2f_circ_mapping: using cached value', 3);
0033 else
0034 
0035     mdl = fix_model(mdl);
0036     switch size(xyzr,1)
0037       case 3; mapping = contained_elems_2d( mdl, xyzr );
0038       case 4; [mapping failed] = contained_elems_3d( mdl, xyzr );
0039       otherwise; error('size of xyzr incorrect');
0040     end
0041 
0042     eidors_obj('set-cache', c_obj, 'circle_mapping', mapping);
0043     eidors_obj('set-cache', c_obj, 'failed_circle_mapping', failed);
0044     eidors_msg('mk_coarse_fine_mapping: setting cached value', 3);
0045 end
0046 
0047 % Mapping depends only on nodes and elems - remove the other stuff
0048 function c_obj = cache_obj(mdl, xyzr)
0049    c_obj = {mdl.nodes, mdl.elems, xyzr};
0050 
0051 % Redirector during test code dev
0052 function mapping = contained_elems_2d( mdl, xyr );
0053 %mapping = contained_elems_2d_new( mdl, xyr );
0054  mapping = contained_elems_2d_old( mdl, xyr );
0055 
0056 function mapping = contained_elems_2d_new( mdl, xyr );
0057    % We fill sparse by columns, (ie adding in CCS storage)
0058    Nc = size(xyr,      2); % Num circs
0059    too_far = elems_too_far( mdl, xyr );
0060 
0061    mapping = sparse( num_elems(mdl) , Nc );
0062    for i=1:Nc
0063      mapping(:,i) = circ_in_elem_2d(mdl, find( ~too_far(:,i)), ...
0064                 xyr(1,i), xyr(2,i), xyr(3,i));
0065    end
0066 
0067 % look only at elements 'look'
0068 function mapping = circ_in_elem_2d( mdl, look, xc, yc, rc);
0069    Nt = elem_dim(mdl) + 1; % nodes per simplex
0070    pirc2 = pi*rc^2;
0071 % start assuming no content
0072    mapping = sparse(num_elems(mdl),1);
0073 % For each element, find nodes inside
0074    els = mdl.elems(look,:);
0075    ndx = reshape(mdl.nodes(els,1) - xc, size(els));
0076    ndy = reshape(mdl.nodes(els,2) - yc, size(els));
0077    n_in = (ndx.^2 + ndy.^2) < rc^2; % node inside
0078 % triangles with 3 nodes inside are all in, stop looking at them
0079    all_n_in = sum(n_in,2) == Nt;
0080    mapping(look(all_n_in)) = 1;
0081    look(all_n_in)= []; n_in(all_n_in,:)= [];
0082 % find distance inside each face
0083    f_in = zeros( length(look), Nt); 
0084    k=1;  for i= look(:)';
0085       faces = mdl.elem2face(i,:);
0086       out   =~mdl.inner_normal(i,:);
0087       f_norm= mdl.normals( faces, :);
0088       f_norm(out,:) = -f_norm(out,:);
0089 
0090       f_ctr = mdl.face_centre( faces,:);
0091       v_ctr = repmat([xc,yc],Nt,1) - f_ctr;
0092       v_ctr = sum(v_ctr .* f_norm,2)/rc; % >1 in, <-1 out,
0093       f_in(k,:) = v_ctr';
0094    k=k+1;end
0095 
0096 % triangles with any sides outside are out
0097    any_s_out= any(f_in<-1,2);
0098    look(any_s_out)= [];
0099    n_in(any_s_out,:) = [];
0100    f_in(any_s_out,:) = [];
0101 
0102 % triangles with 3 sides inside are all in.
0103    all_s_in = sum(f_in>1,2) == Nt;
0104    mapping(look(all_s_in)) = pirc2 / ...
0105             mdl.elem_volume(look(all_s_in));
0106    look(all_s_in)= [];
0107    n_in(all_s_in) = [];
0108    f_in(all_s_in) = [];
0109 
0110 % Now, all triangles should be partially in
0111 % Calculate area chopped out
0112    fin1 = f_in<1;
0113    a_out = zeros(size(fin1));
0114    a_out(fin1) = acos(f_in(fin1));
0115    a_out(fin1) = (a_out(fin1) - cos(a_out(fin1)).*sin(a_out(fin1)))/pi;
0116    
0117    % start with the default. This is accurate if there are
0118    % no contained nodes, otherwise we need to add back
0119    % those fractions
0120    mapping(look) = pirc2 / mdl.elem_volume(look);
0121 
0122    %TODO: rewrite loop to avoid case 0.
0123    k=1; for i= look(:)';
0124       vol = pi*rc^2 / mdl.elem_volume(i);
0125       switch sum(n_in(k,:))
0126          case 0; % already do this
0127    
0128          case 1; 
0129             nd = mdl.elems(k, n_in(k,:));
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    if 0 % for biggish Nc, this is insane
0179        nodes = repmat(mdl.nodes,[1,1,Nc]);
0180        targets = repmat(xyr(1:mdl_dim(mdl),:),[1,1,num_nodes(mdl)]);
0181        targets = shiftdim(targets,2);
0182        dist = nodes - targets;
0183        dist = sqrt(sum(dist.^2,2));
0184        node_target_dist = squeeze(dist);
0185        furthest_elem_node_dist = node_target_dist(mdl.elems,:);
0186        furthest_elem_node_dist = reshape(furthest_elem_node_dist,Ne,Nt,Nc);
0187        [furthest_elem_node_dist, furthest_elem_node]= max(furthest_elem_node_dist,[],2);
0188        furthest_elem_node_dist = squeeze(furthest_elem_node_dist);
0189        furthest_elem_node = squeeze(furthest_elem_node);
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    else
0194        too_far = false(Ne,Nc);
0195       for i = 1:Nc 
0196           targets = repmat(xyr(1:mdl_dim(mdl),i),[1,num_nodes(mdl)])';
0197           dist = mdl.nodes - targets;
0198           dist = sqrt(sum(dist.^2,2));
0199           furthest_elem_node_dist =max(dist(mdl.elems),[],2);
0200           too_far(:,i) = (furthest_elem_node_dist - mdl.max_edge_len) > xyr(Nt,i);
0201       end
0202    end
0203   
0204    
0205   
0206 
0207 function [mapping failed] = contained_elems_3d( mdl, xyr );
0208    Ne = size(mdl.elems,1); % Num elems
0209    Nc = size(xyr,      2); % Num circs
0210    failed(1:Nc) = false;
0211    % We fill sparse by columns, due to CCS storage, this is fairly efficient
0212    mapping = sparse( Ne, Nc );
0213     if 0
0214    % INterpolate
0215    n_interp = 4; % 7-df
0216    m_pts = interp_mesh( mdl, n_interp); 
0217    for i=1:Nc
0218      mapping(:,i) = contained_elem_pts(m_pts, xyr(:,i));
0219    end
0220     else
0221 
0222    % 4. Make a tmp model with only the remaining elems
0223    % 5. Interpolate
0224    % 6. Merge
0225 
0226    too_far = elems_too_far( mdl, xyr );
0227    
0228    tmp = eidors_obj('fwd_model','tmp','nodes',mdl.nodes,'elems',mdl.elems);
0229    %mapping = sparse( Ne, Nc );
0230    n_interp_min = 4;
0231    n_interp_max = 10;
0232    for i=1:Nc
0233        good = ~too_far(:,i);
0234        if ~any(good), continue, end %point outside the mesh
0235        tmp.elems = mdl.elems(good,:);
0236        n_interp = n_interp_min-1;
0237        log_level = eidors_msg('log_level',1);
0238        while(sum(mapping(good,i))==0 && n_interp < n_interp_max-1)
0239            n_interp = n_interp+1;
0240            m_pts = interp_mesh( tmp, n_interp);
0241            mapping(good,i) = contained_elem_pts(m_pts, xyr(:,i));
0242        end
0243        eidors_msg('log_level', log_level);
0244        if (sum(mapping(good,i)) == 0)
0245            failed(i) = true;
0246            eidors_msg(['mk_c2f_circ_mapping: Interpolation failed for point ' num2str(i)]);
0247        end
0248    end
0249     end
0250    
0251    
0252 function frac= contained_elem_pts(m_pts, xyr);
0253 % This is more clear
0254 %    xc = m_pts(:,1,:) - xyr(1);
0255 %    yc = m_pts(:,2,:) - xyr(2);
0256 %    zc = m_pts(:,3,:) - xyr(3);
0257 %    inr= xc.^2 + yc.^2 + zc.^2 < xyr(4)^2;
0258 
0259 % But this is how to stop matlab from wasting memory
0260      inr = (m_pts(:,1,:) - xyr(1)).^2 + ...% xc =
0261            (m_pts(:,2,:) - xyr(2)).^2 + ...% yc =
0262            (m_pts(:,3,:) - xyr(3)).^2;     % zc =
0263      inpts = inr < xyr(4)^2;
0264 
0265 %    frac= mean( inpts ,3);
0266 %    FIXME: Octave doesn't like to mean on logical
0267      frac= mean( int8( inpts ) ,3);
0268      if sum(inpts(:))==0
0269          % TODO: This message is outdated
0270          eidors_msg(['mk_c2f_circ_mapping: Interpolation failed: increase ', ...
0271                          'fwd_model.interp_mesh.n_interp']);
0272      end
0273 
0274 function do_unit_test
0275    %2D example
0276    imdl = mk_common_model('a2c2',16); fmdl=imdl.fwd_model;
0277    xyc = [0,0.27,0.18;0,-0.1,0.03;0,0.1,0.2;0.1,0.37,0.1]';
0278    th=linspace(0,2*pi,20)';
0279    xx=[0*th+1]*xyc(1,:)+sin(th)*xyc(3,:);
0280    yy=[0*th+1]*xyc(2,:)+cos(th)*xyc(3,:);
0281    show_fem(fmdl,[0,0,1]); set(line(xx,yy),'LineWidth',2);
0282 
0283    c2f= mk_c2f_circ_mapping( fmdl, [0;0;0.1] );
0284    t1= all( abs(c2f(1:4)-0.2857)<.001 ) & all( c2f(5:end)==0 );
0285    unit_test_cmp('2D ex 1:',t1,1);
0286 
0287    c2f= mk_c2f_circ_mapping( fmdl, [.0;.05;0.03]); 
0288    t2= abs( c2f(1) - 0.1429) < .001 & all( c2f(2:end)==0 );
0289    unit_test_cmp('2D ex 2:',t2,1);
0290       
0291    %3D example - cylinder
0292    imdl = mk_common_model('a3cr',16); fmdl=imdl.fwd_model;
0293    c2f= mk_c2f_circ_mapping( fmdl, [0;0;0.1]); 
0294    t3= all( abs(c2f(1:4)-0.1714)<.001 ) & all( c2f(5:64)==0 );
0295    unit_test_cmp('3D ex 1:',t2,1);
0296 
0297    %3D example - cylinder
0298    imdl = mk_common_model('a3cr',16); fmdl=imdl.fwd_model;
0299    c2f= mk_c2f_circ_mapping( fmdl, [0;0;0;0.1]); 
0300    unit_test_cmp('3D ex 2a:',c2f(193:196),.0571,.001);
0301    unit_test_cmp('3D ex 2b:',c2f(1:64),0);
0302 
0303    %3D example - cylinder - 2 pts
0304    imdl = mk_common_model('a3cr',16); fmdl=imdl.fwd_model;
0305    c2f= mk_c2f_circ_mapping( fmdl, [0 0;0 0;0 0;0.1 0.2]); 
0306    unit_test_cmp('3D ex 3a:',c2f(193:196),.0571,.001);
0307    unit_test_cmp('3D ex 3b:',c2f(1:64),0);