0001 function [fmdlo]= join_models(fmdl1, fmdl2, tol)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 if ischar(fmdl1) && strcmp(fmdl1,'UNIT_TEST'); do_unit_test; return; end
0017
0018 if nargin==2; tol = 1e-5; end
0019
0020 fmdlo = do_join_models(fmdl1,fmdl2,tol);
0021
0022 function fmdlo = do_join_models(fmdl1,fmdl2,thresh);
0023 n1 = fmdl1.nodes;
0024 n2 = fmdl2.nodes;
0025 Dnodes = 0;
0026 nD = size(n1,2);
0027 for i=1:nD
0028 Dnodes = Dnodes + abs(bsxfun(@minus, n1(:,i), n2(:,i)'));
0029 end
0030
0031 idx = find(any( Dnodes<thresh,1 ));
0032
0033 nn1 = num_nodes(fmdl1);
0034 fmdlo = eidors_obj('fwd_model','joined model');
0035 fmdlo.nodes = [fmdl1.nodes;fmdl2.nodes];
0036 fmdlo.elems = [fmdl1.elems;fmdl2.elems+nn1];
0037 fmdlo.gnd_node = fmdl1.gnd_node;
0038 oidx = 1:num_nodes(fmdlo);
0039 nno = num_nodes(fmdlo);
0040 oidx(idx + nn1)= [];
0041 nidx = zeros(num_nodes(fmdlo),1);
0042 nidx(oidx) = 1:length(oidx);
0043 for i= idx(:)'
0044 ff = find(Dnodes(:,i)<thresh);
0045 if length(ff)~=1; error('Degenerate models with equal nodes'); end
0046 nidx(i+nn1) = ff;
0047 end
0048 fmdlo.nodes = fmdlo.nodes(oidx,:);
0049 fmdlo.elems = reshape(nidx(fmdlo.elems(:)),[],nD+1);
0050 if isfield(fmdl1,'electrode')
0051 fmdlo.electrode = fmdl1.electrode;
0052 else
0053
0054 fmdlo.electrode = struct('nodes',{},'z_contact',{});
0055 end
0056 if isfield(fmdl2,'electrode')
0057 for i=1:length(fmdl2.electrode);
0058 eli = fmdl2.electrode(i);
0059 eli.nodes = nidx(eli.nodes + nn1);
0060 fmdlo.electrode(end+1) = eli;
0061 end
0062 end
0063 fmdlo.boundary = find_boundary(fmdlo.elems);
0064
0065
0066
0067 function do_unit_test
0068
0069 subplot(221); [fmdl1,fmdl2]=do_unit_test_2D_mdls;
0070
0071 fmdlo= join_models(fmdl1, fmdl2);
0072 do_testing('join_models-2Da',fmdlo,fmdl1,fmdl2)
0073 subplot(222); show_fem(fmdlo,[0,1,1]); axis off
0074
0075 fmdl1 = rmfield(fmdl1,'electrode');
0076 fmdlo= join_models(fmdl1, fmdl2);
0077 do_testing('join_models-2Db',fmdlo,fmdl1,fmdl2)
0078
0079 fmdl2 = rmfield(fmdl2,'electrode');
0080 fmdlo= join_models(fmdl1, fmdl2);
0081 do_testing('join_models-2Dc',fmdlo,fmdl1,fmdl2)
0082
0083 subplot(223); [fmdl1,fmdl2]=do_unit_test_3D_mdls; view(0,63);
0084 fmdlo= join_models(fmdl1, fmdl2);
0085 do_testing('join_models-3Da',fmdlo,fmdl1,fmdl2)
0086 subplot(224); show_fem(fmdlo,[0,1,0]); axis off; view(0,63);
0087
0088
0089 function do_testing(txt,fmdlo,fmdl1,fmdl2)
0090 unit_test_cmp([txt,'-01'], ...
0091 size(fmdl1.elems,1) + size(fmdl2.elems,1), ...
0092 size(fmdlo.elems,1) );
0093 unit_test_cmp([txt,'-02'], ...
0094 size(fmdl1.nodes,1) + size(fmdl2.nodes,1) >= ...
0095 size(fmdlo.nodes,1), true );
0096
0097 if ~isfield(fmdl1,'electrode'); fmdl1.electrode = struct([]); end
0098 if ~isfield(fmdl2,'electrode'); fmdl2.electrode = struct([]); end
0099 unit_test_cmp([txt,'-02'], ...
0100 length(fmdl1.electrode) + length(fmdl2.electrode), ...
0101 length(fmdlo.electrode));
0102
0103 function [fmdl1,fmdl2]=do_unit_test_2D_mdls;
0104 imdl = mk_common_model('a2c0',8); fmdl= imdl.fwd_model;
0105 fmdl = crop_model(fmdl,inline('x<-0.4','x','y','z'));
0106 idx = fmdl.nodes(:,1)<-0.25;
0107 fmdl.nodes(idx,1) = -0.25;
0108 fmdl.nodes(:,1) = fmdl.nodes(:,1) + 0.25;
0109 fmdl1 = crop_model(fmdl,inline('x> 1.0','x','y','z'));
0110 fmdl2 = fmdl;
0111 fmdl2.nodes(:,1) = -fmdl2.nodes(:,1);
0112 fmdl2 = crop_model(fmdl2,inline('x+y<-1.28','x','y','z'));
0113 fmdl2 = crop_model(fmdl2,inline('y<-0.95','x','y','z'));
0114 fmdl2 = crop_model(fmdl2,inline('x<-1.20','x','y','z'));
0115 show_fem(fmdl1,[0,1,1]);
0116 hold on; hh=show_fem(fmdl2,[0,1,1]); set(hh,'EdgeColor',[0,0,1]);
0117 hold off; xlim([-1.3,1.3]); axis off
0118
0119 function [fmdl1,fmdl2]=do_unit_test_3D_mdls;
0120 imdl = mk_common_model('n3r2',[16,2]); fmdl= imdl.fwd_model;
0121 fmdl1= crop_model(fmdl,inline('x<-0.55','x','y','z'));
0122 idx = fmdl1.nodes(:,1)<-0.25;
0123 fmdl1.nodes(idx,1) = -0.35;
0124 fmdl1.electrode([19,9])=[];
0125 fmdl1.nodes(:,1) = fmdl1.nodes(:,1) + 0.35;
0126 fmdl2 = fmdl1;
0127 fmdl2.nodes(:,1) = -fmdl2.nodes(:,1);
0128 show_fem(fmdl1,[0,1,0]);
0129 hold on; hh=show_fem(fmdl2,[0,1,0]); set(hh,'EdgeColor',[0,0,1]);
0130 hold off; xlim([-1.3,1.3]); axis off