mk_coarse_fine_mapping

 MK_COARSE_FINE_MAPPING: create a mapping matrix from coarse to fine FEM
 mapping= mk_coarse_fine_mapping( f_mdl, c_mdl );

 Mapping approximates elem_data_fine from elem_data_coase
   elem_data_fine = Mapping*elem_data_coase

 c_mdl is coarse fwd_model
 f_mdl is fine fwd_model

 if the geometry of the fine and coarse models are not
  aligned, then they can be translated and mapped using
    coarse_xyz = M*( fine_xyz - T)
  where
    T= c_mdl.mk_coarse_fine_mapping.f2c_offset (1xN_dims)
    M= c_mdl.mk_coarse_fine_mapping.f2c_project (N_dimsxN_dims)
  by default T= [0,0,0] and M=eye(3)

 if c_mdl is 2D and f_mdl is 3D, then parameter
     c_mdl.mk_coarse_fine_mapping.z_depth
     indicates the +/- z_depth which elements in 2D are
     considered to be extruded in 3D (default inf)

0001 function mapping = mk_coarse_fine_mapping( f_mdl, c_mdl );
0002 % MK_COARSE_FINE_MAPPING: create a mapping matrix from coarse to fine FEM
0003 % mapping= mk_coarse_fine_mapping( f_mdl, c_mdl );
0004 %
0005 % Mapping approximates elem_data_fine from elem_data_coase
0006 %   elem_data_fine = Mapping*elem_data_coase
0007 %
0008 % c_mdl is coarse fwd_model
0009 % f_mdl is fine fwd_model
0010 %
0011 % if the geometry of the fine and coarse models are not
0012 %  aligned, then they can be translated and mapped using
0013 %    coarse_xyz = M*( fine_xyz - T)
0014 %  where
0015 %    T= c_mdl.mk_coarse_fine_mapping.f2c_offset (1xN_dims)
0016 %    M= c_mdl.mk_coarse_fine_mapping.f2c_project (N_dimsxN_dims)
0017 %  by default T= [0,0,0] and M=eye(3)
0018 %
0019 % if c_mdl is 2D and f_mdl is 3D, then parameter
0020 %     c_mdl.mk_coarse_fine_mapping.z_depth
0021 %     indicates the +/- z_depth which elements in 2D are
0022 %     considered to be extruded in 3D (default inf)
0023 
0024 % (C) 2007-2008 Andy Adler. License: GPL version 2 or version 3
0025 % $Id: mk_coarse_fine_mapping.html 2819 2011-09-07 16:43:11Z aadler $
0026 
0027 c_mdl = assign_defaults( c_mdl, f_mdl );
0028 c_obj = cache_obj(c_mdl, f_mdl);
0029 
0030 f_mdl= offset_and_project( f_mdl, c_mdl);
0031 mapping = eidors_obj('get-cache', c_obj, 'coarse_fine_mapping');
0032 if ~isempty(mapping)
0033     eidors_msg('mk_coarse_fine_mapping: using cached value', 3);
0034 else
0035 
0036     z_depth = c_mdl.mk_coarse_fine_mapping.z_depth;
0037 
0038     f_elems = all_contained_elems( f_mdl, c_mdl, z_depth);
0039     mapping = contained_elems_i( f_mdl, c_mdl, f_elems, z_depth);
0040 
0041     if isfield(c_mdl,'coarse2fine')
0042        mapping = mapping*c_mdl.coarse2fine;
0043     end
0044 
0045     eidors_obj('set-cache', c_obj, 'coarse_fine_mapping', mapping);
0046     eidors_msg('mk_coarse_fine_mapping: setting cached value', 3);
0047 end
0048 
0049 % Mapping depends only on nodes and elems - remove the other stuff
0050 function c_obj = cache_obj(c_mdl, f_mdl)
0051    c_obj = {c_mdl.nodes, c_mdl.elems, c_mdl.mk_coarse_fine_mapping, ...
0052             f_mdl.nodes, f_mdl.elems};
0053 
0054 % find all elems of ff_mdl completely contained in cc_mdl
0055 function c_elems = all_contained_elems( fm, cm, z_depth)
0056     [nf,ef]= size(fm.elems);
0057     [nc,ec]= size(cm.elems);
0058     fm_pts = zeros(nf*ef,3);
0059     % shrink pts slightly (s_fac) so they're not on the boundary
0060     % by shrinking, we avoid cases where an element is
0061     % only slighly intersecting another. This is beyond the
0062     % resolution of the next step (interpolation) anyway
0063     s_fac= .9; % .9999;
0064     for dim= 1:ef-1 % ef-1 is dimensions in fm
0065        % fm_pts is local_nodes x elems x xyz
0066        fm_pt= reshape(fm.nodes(fm.elems,dim),nf,ef);
0067        fm_ctr= mean(fm_pt,2)*ones(1,ef);
0068        fm_pt = s_fac*(fm_pt-fm_ctr) + fm_ctr;
0069        fm_pts(:,dim) = fm_pt(:);
0070     end
0071 
0072     tsn= search_fm_pts_in_cm(cm, fm_pts, z_depth);
0073     tsn= reshape( tsn, [], ef);
0074     % if all points are outside (NaN) then c_elems = -1
0075     % if all points are in one elem   then c_elems = elem #
0076     % if all points are in diff elems then c_elems = 0
0077     c_elems= all(diff(tsn,1,2)==0,2) .* tsn(:,1);
0078     c_elems(all(tsn==-1,2))= -1; % all points outside
0079 
0080 
0081 % tsn = vector of length size(fm_pts,1)
0082 % tsn(i) = elem in cm which contains point
0083 % tsn(i) = -1 if point is outside cm (and z_depth, if appropriate)
0084 function tsn= search_fm_pts_in_cm(cm, fm_pts, z_depth);
0085     dc= size(cm.elems,2)-1;  %coarse dim
0086     df= size(fm_pts,2); %fine dim
0087 
0088     tsn= -ones(size(fm_pts,1),1);
0089     not_oor= (tsn==-1); % logical 1
0090 
0091     if dc==2  %corse model is 2D
0092 
0093        if df==3
0094        % look for f_mdl z not out of range
0095           not_oor= not_oor &  any( abs(fm_pts(:,3) ) <= z_depth , 2);
0096        end
0097        dims=1:2;
0098 
0099     elseif dc==3 %coarse model is 3D
0100 
0101        dims=1:3; 
0102 
0103     else
0104        error('coarse model must be 2 or 3D');
0105     end
0106 
0107     tsn(not_oor)= tsearchn(cm.nodes(:,dims), cm.elems, fm_pts(not_oor,dims));
0108     tsn(isnan(tsn))= -1;
0109 
0110 % find fraction of elem contained in cm
0111 % idx is the index into which the elem is contained
0112 % if idx >= 1, the element is completely with in that coarse elem
0113 % if idx == 0, the element is crosses several coarse elems
0114 % if idx ==-1, the element is outside the coarse model
0115 function c_elems = contained_elems_i( fm, cm, idx, z_depth)
0116    [nc,dc]= size(cm.elems);
0117    [nf,df]= size(fm.elems);
0118 
0119    fidx= find(idx==0);
0120    l_fidx= length(fidx);
0121 
0122    c_e_i= []; c_e_j=[]; c_e_v=[];
0123 
0124    % lower density interpolation in higher dimentions, since
0125    % the added dimensions will give extra interpolation points.
0126    n_interp = 7-df;
0127    interp_mdl.nodes= fm.nodes;
0128    interp_mdl.elems= fm.elems(fidx,:);
0129 
0130    fm_pts = interp_mesh( interp_mdl, n_interp);
0131    l_interp = size(fm_pts,3);
0132 
0133    fm_pts = permute( fm_pts, [3,1,2]);
0134    fm_pts = reshape(fm_pts, [], df-1);
0135 
0136    tsn= search_fm_pts_in_cm(cm, fm_pts, z_depth);
0137 
0138    tsn_idx= ones(l_interp,1)*fidx(:)';
0139    tsn_idx= tsn_idx(:);
0140    % find and isolate outside elements
0141    outside_idx= tsn==-1;
0142    tsn(outside_idx) = [];
0143    tsn_idx(outside_idx) = [];
0144    
0145    in_idx= find(idx<=0);
0146    ridx= 1:nf; ridx(in_idx)=[];
0147    idx(in_idx)=[];
0148 
0149    % first term is contribution from f_elems in one c_elem
0150    % next term is contribution from f_elems in many c_elems, weighted
0151    c_elems = sparse(ridx,idx,1,nf,nc) +  ...
0152              sparse(tsn_idx,tsn,1,nf,nc)/l_interp;
0153       
0154 % Do 3D interpolation of region xyzmin= [x,y,z] to xyzmax
0155 %  with n_interp points in the minimum direction
0156 function xyz = interpxyz( xyzmin, xyzmax, n_interp);
0157     xyzdelta= xyzmax - xyzmin;
0158     xyz_interp = 1 + floor(n_interp * xyzdelta / max(xyzdelta) );
0159     xspace = linspace(xyzmin(1), xyzmax(1), xyz_interp(1) );
0160     yspace = linspace(xyzmin(2), xyzmax(2), xyz_interp(2) );
0161     zspace = linspace(xyzmin(3), xyzmax(3), xyz_interp(3) );
0162     [xx3,yy3,zz3] = ndgrid( xspace, yspace, zspace );
0163     xyz= [xx3(:), yy3(:), zz3(:)];
0164 
0165 % Offset and project f_mdl as required
0166 function f_mdl= offset_and_project( f_mdl, c_mdl);
0167     [fn,fd]= size(f_mdl.nodes);
0168     T= c_mdl.mk_coarse_fine_mapping.f2c_offset;
0169     M= c_mdl.mk_coarse_fine_mapping.f2c_project;
0170 
0171     f_mdl.nodes= ( f_mdl.nodes - ones(fn,1)*T )*M;
0172 
0173 function c_mdl = assign_defaults( c_mdl, f_mdl );
0174     [fn,fd]= size(f_mdl.nodes);
0175     try    c_mdl.mk_coarse_fine_mapping.f2c_offset; % test exist
0176     catch  c_mdl.mk_coarse_fine_mapping.f2c_offset= zeros(1,fd);
0177     end
0178     try    c_mdl.mk_coarse_fine_mapping.f2c_project;
0179     catch  c_mdl.mk_coarse_fine_mapping.f2c_project= speye(fd);
0180     end
0181     try    c_mdl.mk_coarse_fine_mapping.z_depth;
0182     catch  c_mdl.mk_coarse_fine_mapping.z_depth= inf;
0183     end