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.html 2819 2011-09-07 16:43:11Z aadler $
0033 
0034 if isstr(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test;return; end
0035 
0036 if nargin<2; n_interp=0; end
0037 try n_interp = mdl.interp_mesh.n_interp; end % Override if provided
0038 
0039 
0040 % cashing
0041    
0042    c_obj = {mdl.elems, mdl.nodes, n_interp};
0043    mdl_pts = eidors_obj('get-cache', c_obj, 'interpolation');
0044    if ~isempty(mdl_pts)
0045        return
0046    end
0047 
0048 
0049 % Get element nodes, and reshape
0050 % need to be of size n_dims_1 x (n_elems*n_dims) for reshape
0051 el_nodes= mdl.nodes(mdl.elems',:);
0052 el_nodes= reshape(el_nodes, elem_dim(mdl)+1, []);
0053 
0054 % Get interpolation matrix
0055 interp= triangle_interpolation( n_interp, elem_dim(mdl) );
0056 l_interp = size(interp,1);
0057 
0058 mdl_pts = interp*el_nodes;
0059 mdl_pts = reshape(mdl_pts, l_interp, num_elems(mdl), mdl_dim(mdl));
0060 
0061 mdl_pts = permute(mdl_pts, [2,3,1]);
0062 
0063 % caching
0064 eidors_cache('boost_priority', -2); % low priority
0065 c_obj = {mdl.elems, mdl.nodes, n_interp};
0066 eidors_obj('set-cache', c_obj, 'interpolation', mdl_pts);
0067 eidors_cache('boost_priority', +2); % restore priority
0068 
0069 
0070 % interpolate over a triangle with n_interp points
0071 % generate a set of points to fairly cover the triangle
0072 % dim_coarse is dimensions + 1 of coarse model
0073 function interp= triangle_interpolation(n_interp, el_dim)
0074     interp= zeros(0,el_dim+1);
0075 
0076     switch el_dim
0077      case 1;
0078        for i=0:n_interp
0079              interp= [interp;i,n_interp-i];
0080        end
0081      case 2;
0082        for i=0:n_interp
0083           for j=0:n_interp-i
0084              interp= [interp;i,j,n_interp-i-j];
0085           end
0086        end
0087      case 3;
0088        for i=0:n_interp
0089           for j=0:n_interp-i
0090              for k=0:n_interp-i-j
0091                 interp= [interp;i,j,k,n_interp-i-j-k];
0092              end
0093           end
0094        end
0095      otherwise;
0096        error('Element is not 1D (line), 2D (triangle) or 3D (tetrahedron)');
0097     end
0098 
0099     interp= (interp + 1/(el_dim+1) )/(n_interp+1);
0100 
0101 
0102 function do_unit_test
0103     mdl.type='fwd_model';mdl.name='jnk';
0104     mdl.nodes= [4,6,8,4]';
0105     mdl.elems=[1,2;2,4;2,3;4,4];
0106     pp=interp_mesh(mdl,0);
0107     unit_test_cmp('1D/1D (#1): ',pp,[5;5;7;4]);
0108 
0109     mdl.nodes= [4,6,8,4;2,2,2,6]';
0110     pp=interp_mesh(mdl,0);
0111     unit_test_cmp('1D/2D (#1): ',pp,[5 2;5 4;7 2;4 6]);
0112 
0113     mdl.nodes= 3*[4,6,8,4,6,8;2,2,2,5,5,5]';
0114     mdl.elems=[1,2,4;2,4,5;2,3,5;3,5,6];
0115     pp=interp_mesh(mdl,0);
0116     unit_test_cmp('2D/2D (#1): ',pp,[14 9;16 12; 20 9;22 12]);
0117 
0118     pp=interp_mesh(mdl,3);
0119     unit_test_cmp('2D/2D (#2): ',pp(:,:,6),[14 9;16 12; 20 9;22 12]);
0120 
0121     mdl.nodes = [mdl.nodes, 0*mdl.nodes(:,1)+3];
0122     pp=interp_mesh(mdl,0);
0123     unit_test_cmp('2D/3D (#1): ',pp,[14 9 3;16 12 3; 20 9 3;22 12 3]);
0124 
0125     mdl = mk_common_model('n3r2',16); mdl= mdl.fwd_model;
0126     pp=interp_mesh(mdl,0);
0127     unit_test_cmp('3D/3D (#1): ',pp(1,:),[0.920196320100808   0.048772580504032   0.5],1e-14);
0128 
0129     pp=interp_mesh(mdl,4);
0130     unit_test_cmp('3D/3D (#2a):',pp(1,:,21),[0.920196320100808   0.048772580504032   0.5],1e-14);
0131