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 5553 2017-06-18 16:35:36Z aadler $
0022 
0023 if ischar(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0024 
0025 copt.cache_obj = cache_obj(mdl, xyzr);
0026 copt.fstr = 'mk_c2f_circ_mapping';
0027 [mapping, failed] = eidors_cache(@circ_mapping,{mdl,xyzr},copt);
0028 
0029 function [mapping, failed] = circ_mapping(mdl,xyzr,copt)
0030 
0031     failed = false;% not all subfunctions set this
0032     mdl = fix_model(mdl);
0033     switch size(xyzr,1)
0034       case 3; % use for 2D or for cylinder mapping in 3D
0035          mapping = contained_elems_2d( mdl, xyzr );
0036          if mdl_dim(mdl) == 2;
0037             correctmap = pi*xyzr(end,:).^2;
0038          else
0039             % No analytic way to calculate correct value (for non extruded models)
0040             correctmap = (get_elem_volume(mdl)'*mapping);
0041          end
0042       case 4; % use for 3D and spherical maps
0043          [mapping failed] = contained_elems_3d( mdl, xyzr );
0044          correctmap = 4/3*pi*xyzr(end,:).^3;
0045       otherwise; error('size of xyzr incorrect');
0046     end
0047 
0048     % Correct
0049     vol = get_elem_volume(mdl,-2)'; %-2 => don't use c2f
0050 
0051     mapping = bsxfun(@times, mapping, correctmap./(vol*mapping));
0052     
0053 % Mapping depends only on nodes and elems - remove the other stuff
0054 function c_obj = cache_obj(mdl, xyzr)
0055    c_obj = {mdl.nodes, mdl.elems, xyzr};
0056 
0057 % Redirector during test code dev
0058 function mapping = contained_elems_2d( mdl, xyr );
0059 %mapping = contained_elems_2d_new( mdl, xyr );
0060  mapping = contained_elems_2d_old( mdl, xyr );
0061 
0062 function mapping = contained_elems_2d_new( mdl, xyr );
0063    % We fill sparse by columns, (ie adding in CCS storage)
0064    Nc = size(xyr,      2); % Num circs
0065    too_far = elems_too_far( mdl, xyr );
0066 
0067    mapping = sparse( num_elems(mdl) , Nc );
0068    for i=1:Nc
0069      mapping(:,i) = circ_in_elem_2d(mdl, find( ~too_far(:,i)), ...
0070                 xyr(1,i), xyr(2,i), xyr(3,i));
0071    end
0072 
0073 % look only at elements 'look'
0074 function mapping = circ_in_elem_2d( mdl, look, xc, yc, rc);
0075    Nt = elem_dim(mdl) + 1; % nodes per simplex
0076    pirc2 = pi*rc^2;
0077 % start assuming no content
0078    mapping = sparse(num_elems(mdl),1);
0079 % For each element, find nodes inside
0080    els = mdl.elems(look,:);
0081    ndx = reshape(mdl.nodes(els,1) - xc, size(els));
0082    ndy = reshape(mdl.nodes(els,2) - yc, size(els));
0083    n_in = (ndx.^2 + ndy.^2) < rc^2; % node inside
0084 % triangles with 3 nodes inside are all in, stop looking at them
0085    all_n_in = sum(n_in,2) == Nt;
0086    mapping(look(all_n_in)) = 1;
0087    look(all_n_in)= []; n_in(all_n_in,:)= [];
0088 % find distance inside each face
0089    f_in = zeros( length(look), Nt); 
0090    k=1;  for i= look(:)';
0091       faces = mdl.elem2face(i,:);
0092       out   =~mdl.inner_normal(i,:);
0093       f_norm= mdl.normals( faces, :);
0094       f_norm(out,:) = -f_norm(out,:);
0095 
0096       f_ctr = mdl.face_centre( faces,:);
0097       v_ctr = repmat([xc,yc],Nt,1) - f_ctr;
0098       v_ctr = sum(v_ctr .* f_norm,2)/rc; % >1 in, <-1 out,
0099       f_in(k,:) = v_ctr';
0100    k=k+1;end
0101 
0102 % triangles with any sides outside are out
0103    any_s_out= any(f_in<-1,2);
0104    look(any_s_out)= [];
0105    n_in(any_s_out,:) = [];
0106    f_in(any_s_out,:) = [];
0107 
0108 % triangles with 3 sides inside are all in.
0109    all_s_in = sum(f_in>1,2) == Nt;
0110    mapping(look(all_s_in)) = pirc2 / ...
0111             mdl.elem_volume(look(all_s_in));
0112    look(all_s_in)= [];
0113    n_in(all_s_in) = [];
0114    f_in(all_s_in) = [];
0115 
0116 % Now, all triangles should be partially in
0117 % Calculate area chopped out
0118    fin1 = f_in<1;
0119    a_out = zeros(size(fin1));
0120    a_out(fin1) = acos(f_in(fin1));
0121    a_out(fin1) = (a_out(fin1) - cos(a_out(fin1)).*sin(a_out(fin1)))/pi;
0122    
0123    % start with the default. This is accurate if there are
0124    % no contained nodes, otherwise we need to add back
0125    % those fractions
0126    mapping(look) = pirc2 / mdl.elem_volume(look);
0127 
0128    %TODO: rewrite loop to avoid case 0.
0129    k=1; for i= look(:)';
0130       vol = pi*rc^2 / mdl.elem_volume(i);
0131       switch sum(n_in(k,:))
0132          case 0; % already do this
0133    
0134          case 1; 
0135             nd = mdl.elems(k, n_in(k,:));
0136             vol = vol + pi_slice(p1,p2,[xc,yc],mdl.nodes(nd,:),rc);
0137 
0138          case 2; 
0139             nd = mdl.elems(k, n_in(k,:));
0140             vol = vol ...
0141              + pi_slice(p1,p2,[xc,yc],mdl.nodes(nd(1),:),rc) ...
0142              + pi_slice(p1,p2,[xc,yc],mdl.nodes(nd(2),:),rc);
0143 
0144          otherwise; error('cant get here'); 
0145       end
0146    k=k+1; end
0147    
0148 
0149 % Calculate the area of a slice
0150 function a = pi_slice(p1,p2,c,p,r)
0151   a_p12 = 0.5*abs(det([1,p1;1,p2;1,p]));
0152 
0153   a_c12 = 0.5*abs(det([1,p1;1,p2;1,c]));
0154   np1c  = p1-c; np1c = np1c / norm(np1c);
0155   np2c  = p2-c; np2c = np2c / norm(np2c);
0156   ang   = acos( dot(np1c,np2c) );
0157   area  = ang*r^2/2 - a_c12 + a_p12;
0158      
0159 
0160 function mapping = contained_elems_2d_old( mdl, xyr );
0161    Ne = size(mdl.elems,1); % Num elems
0162    Nc = size(xyr,      2); % Num circs
0163    % We fill sparse by columns, due to CCS storage, this is fairly efficient
0164    mapping = sparse( Ne, Nc );
0165 
0166    % Interpolate
0167    n_interp =  7-size(mdl.nodes,2);
0168    m_pts = interp_mesh( mdl, n_interp); 
0169    for i=1:Nc
0170      xc = m_pts(:,1,:) - xyr(1,i);
0171      yc = m_pts(:,2,:) - xyr(2,i);
0172      inr= xc.^2 + yc.^2 < xyr(3,i)^2;
0173      frac= mean(inr,3);
0174      mapping(:,i) = frac;
0175    end
0176 
0177    % 1. Get furthest node in each element
0178    % 2. Get the max edge length of each elem
0179    % 3. Find elems that are too far
0180 function too_far = elems_too_far( mdl, xyr );
0181    Ne = num_elems(mdl);
0182    Nc = size(xyr, 2); % Num circs
0183    Nt = elem_dim(mdl) + 1; % Elements per simplex
0184    if 0 % for biggish Nc, this is insane
0185        nodes = repmat(mdl.nodes,[1,1,Nc]);
0186        targets = repmat(xyr(1:mdl_dim(mdl),:),[1,1,num_nodes(mdl)]);
0187        targets = shiftdim(targets,2);
0188        dist = nodes - targets;
0189        dist = sqrt(sum(dist.^2,2));
0190        node_target_dist = squeeze(dist);
0191        furthest_elem_node_dist = node_target_dist(mdl.elems,:);
0192        furthest_elem_node_dist = reshape(furthest_elem_node_dist,Ne,Nt,Nc);
0193        [furthest_elem_node_dist, furthest_elem_node]= max(furthest_elem_node_dist,[],2);
0194        furthest_elem_node_dist = squeeze(furthest_elem_node_dist);
0195        furthest_elem_node = squeeze(furthest_elem_node);
0196        max_edge_len =  repmat(mdl.max_edge_len,1,Nc);
0197        radius = ones(Ne,1)*xyr(Nt,:);
0198        too_far = (furthest_elem_node_dist - max_edge_len) > radius;
0199    else
0200        too_far = false(Ne,Nc);
0201        progress_msg('mk_c2f_circ_mapping: prepare models',0,Nc);
0202        for i = 1:Nc
0203           progress_msg(i,Nc);
0204           targets = repmat(xyr(1:mdl_dim(mdl),i),[1,num_nodes(mdl)])';
0205           dist = mdl.nodes - targets;
0206           dist = sqrt(sum(dist.^2,2));
0207           furthest_elem_node_dist =max(dist(mdl.elems),[],2);
0208           too_far(:,i) = (furthest_elem_node_dist - mdl.max_edge_len) > xyr(Nt,i);
0209        end
0210        progress_msg(Inf);
0211    end
0212   
0213    
0214   
0215 
0216 function [mapping failed] = contained_elems_3d( mdl, xyr );
0217    Ne = size(mdl.elems,1); % Num elems
0218    Nc = size(xyr,      2); % Num circs
0219    failed(1:Nc) = false;
0220    % We fill sparse by columns, due to CCS storage, this is fairly efficient
0221    mapping = sparse( Ne, Nc );
0222 
0223        % 4. Make a tmp model with only the remaining elems
0224        % 5. Interpolate
0225        % 6. Merge
0226        
0227        too_far = elems_too_far( mdl, xyr );
0228        
0229        tmp = eidors_obj('fwd_model','tmp','nodes',mdl.nodes,'elems',mdl.elems);
0230        %mapping = sparse( Ne, Nc );
0231        try   
0232            n_interp_min = mdl.interp_mesh.n_points;
0233        catch
0234            n_interp_min = 4;
0235        end
0236        n_interp_max = 10;
0237 
0238    if 0
0239        % INterpolate
0240        n_interp = 4; % 7-df
0241        m_pts = interp_mesh( mdl, n_interp); 
0242        for i=1:Nc
0243          mapping(:,i) = contained_elem_pts(m_pts, xyr(:,i));
0244        end
0245    else
0246 
0247        progress_msg('mk_c2f_circ_mapping: calculate mapping',0,Nc);
0248        for i=1:Nc
0249            progress_msg(i,Nc);
0250            good = ~too_far(:,i);
0251            if ~any(good), continue, end %point outside the mesh
0252            tmp.elems = mdl.elems(good,:);
0253            n_interp = n_interp_min-1;
0254            log_level = eidors_msg('log_level',1);
0255            while(sum(mapping(good,i))==0 && n_interp < n_interp_max-1)
0256                n_interp = n_interp+1;
0257                m_pts = interp_mesh( tmp, n_interp);
0258                mapping(good,i) = contained_elem_pts(m_pts, xyr(:,i));
0259            end
0260            eidors_msg('log_level', log_level);
0261            if (sum(mapping(good,i)) == 0)
0262                failed(i) = true;
0263                eidors_msg(['mk_c2f_circ_mapping: Interpolation failed for point ' num2str(i)]);
0264            end
0265        end
0266        progress_msg(Inf);
0267    end
0268    
0269    
0270 function frac= contained_elem_pts(m_pts, xyr);
0271 % This is more clear
0272 %    xc = m_pts(:,1,:) - xyr(1);
0273 %    yc = m_pts(:,2,:) - xyr(2);
0274 %    zc = m_pts(:,3,:) - xyr(3);
0275 %    inr= xc.^2 + yc.^2 + zc.^2 < xyr(4)^2;
0276 
0277 % But this is how to stop matlab from wasting memory
0278      inr = (m_pts(:,1,:) - xyr(1)).^2 + ...% xc =
0279            (m_pts(:,2,:) - xyr(2)).^2 + ...% yc =
0280            (m_pts(:,3,:) - xyr(3)).^2;     % zc =
0281      inpts = inr < xyr(4)^2;
0282 
0283      frac= mean( inpts ,3);
0284 %    FIXME: Octave doesn't like to mean on logical
0285 %    frac= mean( int8( inpts ) ,3);
0286      if sum(inpts(:))==0
0287          % TODO: This message is outdated
0288          eidors_msg(['mk_c2f_circ_mapping: Interpolation failed: increase ', ...
0289                          'fwd_model.interp_mesh.n_interp']);
0290      end
0291 
0292 function do_unit_test
0293    fmdl = ng_mk_cyl_models([0,1,.1],[16,1],.03);
0294    vol = get_elem_volume(fmdl)';
0295    xyr = ones(3,1)*linspace(-.5,.5,7);
0296    rr = .1; VV = pi*rr^2; xyr(3,:) = rr;
0297    [c2f,fail] = mk_c2f_circ_mapping(fmdl,xyr);
0298    unit_test_cmp('2D #1.1:',sum(fail),0);
0299    unit_test_cmp('2D #1.2(r=.1):',vol*c2f/VV,1,1e-2);
0300 
0301    fmdl = ng_mk_cyl_models([2,1,.1],[16,1],.03);
0302    fmdl.nodes(:,3) = fmdl.nodes(:,3) - 1;
0303    vol = get_elem_volume(fmdl)';
0304    xyzr = ones(4,1)*linspace(-.5,.5,7);
0305 
0306    rr = .1; VV = pi*4/3*rr^3; xyzr(4,:) = rr;
0307    [c2f,fail] = mk_c2f_circ_mapping(fmdl,xyzr);
0308    unit_test_cmp('3D #1.1:',sum(fail),0);
0309    unit_test_cmp('3D #1.2(r=.1):',vol*c2f/VV,1,1e-2);
0310 
0311    rr = .01; VV = pi*4/3*rr^3; xyzr(4,:) = rr;
0312    [c2f,fail] = mk_c2f_circ_mapping(fmdl,xyzr);
0313    unit_test_cmp('3D #1.1(r=.01):',sum(fail),0);
0314    unit_test_cmp('3D #1.2:',vol*c2f/VV,1,1e-2);
0315 
0316    %2D example
0317    imdl = mk_common_model('a2c2',16); fmdl=imdl.fwd_model;
0318    xyc = [0,0.27,0.18;0,-0.1,0.03;0,0.1,0.2;0.1,0.37,0.1]';
0319    th=linspace(0,2*pi,20)';
0320    xx=[0*th+1]*xyc(1,:)+sin(th)*xyc(3,:);
0321    yy=[0*th+1]*xyc(2,:)+cos(th)*xyc(3,:);
0322    show_fem(fmdl,[0,0,1]); set(line(xx,yy),'LineWidth',2);
0323 
0324    % split over four elements
0325    rr= 0.1;c2f= mk_c2f_circ_mapping( fmdl, [0;0;rr] );
0326    tt= zeros(size(c2f)); tt(1:4) = pi*rr^2/4; tt= tt./get_elem_volume(fmdl);
0327    unit_test_cmp('2D ex 1:',c2f,tt,1e-10);
0328 
0329    % all in element #1
0330    rr= 0.03;c2f= mk_c2f_circ_mapping( fmdl, [.0;.05;rr]); 
0331    tt= zeros(size(c2f)); tt(1) = pi*rr^2; tt= tt./get_elem_volume(fmdl);
0332    unit_test_cmp('2D ex 2:',c2f,tt,1e-10);
0333       
0334    %3D example - cylinder
0335    imdl = mk_common_model('a3cr',16); fmdl=imdl.fwd_model;
0336    fmdl.nodes = 1.1*fmdl.nodes;
0337    rr=0.1;c2f= mk_c2f_circ_mapping( fmdl, [0;0;rr]); 
0338    V = pi*rr^2*(max(fmdl.nodes(:,3)) - min(fmdl.nodes(:,3)));
0339    unit_test_cmp('3D ex 1 (cylinder):',get_elem_volume(fmdl)'*c2f,V,1e-2);
0340 
0341    %3D example - sphere
0342    imdl = mk_common_model('a3cr',16); fmdl=imdl.fwd_model;
0343    rr=0.05;c2f= mk_c2f_circ_mapping( fmdl, [0;0;0;rr]); 
0344    tt = 4/3*pi*rr^3/24./get_elem_volume(fmdl);
0345    unit_test_cmp('3D ex 2a:',c2f(193:196),tt(193:196),1e-10);
0346    unit_test_cmp('3D ex 2b:',c2f(1:64),0);
0347 
0348    imdl = mk_common_model('a3cr',16); fmdl=imdl.fwd_model;
0349    rr=0.05;c2f= mk_c2f_circ_mapping( fmdl, [0 0;0 0;0 0;rr,rr]); 
0350    unit_test_cmp('3D ex 3a:',c2f(193:196,:),tt(193:196)*[1,1],1e-10);
0351    unit_test_cmp('3D ex 3b:',c2f(1:64,:),0);