interp_mesh

 INTERP_MESH: calculate interpolation points onto mdl elements
    mdl_pts = interp_mesh( fwd_model, n_interp)
 INPUT:
    fwd_model: fwd_model structure
    n_interp:  order of interpolation
      n_interp = 0 - output elem centres (default)
      n_interp >=1 - output multiple points per elem
           in 2D: (n_int+1)*(n_int+2)/2 points per elem
           in 3D: (n_int+1)*(n_int+2)*(n_int+3)/6 points per elem
   n_interp may be specified as:
      fwd_model.interp_mesh.n_interp (This overrides the above)

 OUTPUT:
    mdl_pts = N_elems x N_dims x N_points
   example: for mdl_pts= interp_mesh( mdl, 0);
           mdl_pts(i,:) is centre of element #i      
   example: for mdl_pts= interp_mesh( mdl, 2);
           mdl_pts(i,:,:) is 1 x [x,y,{z}] x n_points to interpolate
 
 EXAMPLE:
   mdl.nodes= [4,6,8,4,6,8;2,2,2,5,5,5]';
   mdl.elems=[1,2,4;2,4,5;2,3,5;3,5,6];
   mdl.type='fwd_model';mdl.name='jnk';
   pp=interp_mesh(mdl,4);
   show_fem(mdl); hold on ;
      for i=1:size(pp,3); plot(pp(:,1,i),pp(:,2,i),'*'); end ;
   hold off

0001 function mdl_pts = interp_mesh( mdl, n_interp)
0002 % INTERP_MESH: calculate interpolation points onto mdl elements
0003 %    mdl_pts = interp_mesh( fwd_model, n_interp)
0004 % INPUT:
0005 %    fwd_model: fwd_model structure
0006 %    n_interp:  order of interpolation
0007 %      n_interp = 0 - output elem centres (default)
0008 %      n_interp >=1 - output multiple points per elem
0009 %           in 2D: (n_int+1)*(n_int+2)/2 points per elem
0010 %           in 3D: (n_int+1)*(n_int+2)*(n_int+3)/6 points per elem
0011 %   n_interp may be specified as:
0012 %      fwd_model.interp_mesh.n_interp (This overrides the above)
0013 %
0014 % OUTPUT:
0015 %    mdl_pts = N_elems x N_dims x N_points
0016 %   example: for mdl_pts= interp_mesh( mdl, 0);
0017 %           mdl_pts(i,:) is centre of element #i
0018 %   example: for mdl_pts= interp_mesh( mdl, 2);
0019 %           mdl_pts(i,:,:) is 1 x [x,y,{z}] x n_points to interpolate
0020 %
0021 % EXAMPLE:
0022 %   mdl.nodes= [4,6,8,4,6,8;2,2,2,5,5,5]';
0023 %   mdl.elems=[1,2,4;2,4,5;2,3,5;3,5,6];
0024 %   mdl.type='fwd_model';mdl.name='jnk';
0025 %   pp=interp_mesh(mdl,4);
0026 %   show_fem(mdl); hold on ;
0027 %      for i=1:size(pp,3); plot(pp(:,1,i),pp(:,2,i),'*'); end ;
0028 %   hold off
0029 %
0030 
0031 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0032 % $Id: interp_mesh.m 5866 2018-12-05 17:09:40Z aadler $
0033 
0034 if ischar(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test;return; end
0035 
0036 if nargin<2; n_interp=0; end
0037 if strcmp(mdl.type,'image') || strcmp(mdl.type,'inv_model')
0038    mdl = mdl.fwd_model;
0039 end
0040 
0041 
0042 try n_interp = mdl.interp_mesh.n_interp; end % Override if provided
0043 
0044 
0045 copt.cache_obj = {mdl.elems, mdl.nodes, n_interp};
0046 copt.fstr = 'interp_mesh';
0047 copt.boost_priority = -2;
0048 
0049 mdl_pts = eidors_cache(@do_interpolation,{mdl, n_interp},copt);
0050 
0051 
0052 function mdl_pts = do_interpolation(mdl, n_interp)
0053    % Get element nodes, and reshape
0054    % need to be of size n_dims_1 x (n_elems*n_dims) for reshape
0055    el_nodes= mdl.nodes(mdl.elems',:);
0056    el_nodes= reshape(el_nodes, elem_dim(mdl)+1, []);
0057 
0058    % Get interpolation matrix
0059    interp= triangle_interpolation( n_interp, elem_dim(mdl) );
0060    l_interp = size(interp,1);
0061 
0062    mdl_pts = interp*el_nodes;
0063    mdl_pts = reshape(mdl_pts, l_interp, num_elems(mdl), mdl_dim(mdl));
0064 
0065    mdl_pts = permute(mdl_pts, [2,3,1]);
0066    
0067    
0068    
0069 % interpolate over a triangle with n_interp points
0070 % generate a set of points to fairly cover the triangle
0071 % dim_coarse is dimensions + 1 of coarse model
0072 function interp= triangle_interpolation(n_interp, el_dim)
0073     interp= zeros(0,el_dim+1);
0074 
0075     switch el_dim
0076      case 1;
0077        for i=0:n_interp
0078              interp= [interp;i,n_interp-i];
0079        end
0080      case 2;
0081        for i=0:n_interp
0082           for j=0:n_interp-i
0083              interp= [interp;i,j,n_interp-i-j];
0084           end
0085        end
0086      case 3;
0087        for i=0:n_interp
0088           for j=0:n_interp-i
0089              for k=0:n_interp-i-j
0090                 interp= [interp;i,j,k,n_interp-i-j-k];
0091              end
0092           end
0093        end
0094      otherwise;
0095        error('Element is not 1D (line), 2D (triangle) or 3D (tetrahedron)');
0096     end
0097 
0098     interp= (interp + 1/(el_dim+1) )/(n_interp+1);
0099 
0100 
0101 function do_unit_test
0102     mdl.type='fwd_model';mdl.name='jnk';
0103     mdl.nodes= [4,6,8,4]';
0104     mdl.elems=[1,2;2,4;2,3;4,4];
0105     pp=interp_mesh(mdl,0);
0106     unit_test_cmp('1D/1D (#1): ',pp,[5;5;7;4]);
0107 
0108     mdl.nodes= [4,6,8,4;2,2,2,6]';
0109     pp=interp_mesh(mdl,0);
0110     unit_test_cmp('1D/2D (#1): ',pp,[5 2;5 4;7 2;4 6]);
0111 
0112     mdl.nodes= 3*[4,6,8,4,6,8;2,2,2,5,5,5]';
0113     mdl.elems=[1,2,4;2,4,5;2,3,5;3,5,6];
0114     pp=interp_mesh(mdl,0);
0115     unit_test_cmp('2D/2D (#1): ',pp,[14 9;16 12; 20 9;22 12]);
0116 
0117     pp=interp_mesh(mdl,3);
0118     unit_test_cmp('2D/2D (#2): ',pp(:,:,6),[14 9;16 12; 20 9;22 12]);
0119 
0120     mdl.nodes = [mdl.nodes, 0*mdl.nodes(:,1)+3];
0121     pp=interp_mesh(mdl,0);
0122     unit_test_cmp('2D/3D (#1): ',pp,[14 9 3;16 12 3; 20 9 3;22 12 3]);
0123 
0124     mdl = mk_common_model('n3r2',[16,2]); mdl= mdl.fwd_model;
0125     pp=interp_mesh(mdl,0);
0126     unit_test_cmp('3D/3D (#1): ',pp(1,:),[0.920196320100808   0.048772580504032   0.5],1e-14);
0127 
0128     pp=interp_mesh(mdl,4);
0129     unit_test_cmp('3D/3D (#2a):',pp(1,:,21),[0.920196320100808   0.048772580504032   0.5],1e-14);
0130