interp_mesh

PURPOSE ^

INTERP_MESH: calculate interpolation points onto mdl elements

SYNOPSIS ^

function mdl_pts = interp_mesh( mdl, n_interp)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005