0001 function mdl_pts = interp_mesh( mdl, n_interp)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
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
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
0049
0050 el_nodes= mdl.nodes(mdl.elems',:);
0051 el_nodes= reshape(el_nodes, elem_dim(mdl)+1, []);
0052
0053
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
0065
0066
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