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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
