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 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
0038
0039
0040
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
0050
0051 el_nodes= mdl.nodes(mdl.elems',:);
0052 el_nodes= reshape(el_nodes, elem_dim(mdl)+1, []);
0053
0054
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
0064 eidors_cache('boost_priority', -2);
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);
0068
0069
0070
0071
0072
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