Description of mk_c2f_circ

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
