1 function mesh=Ujoinmeshes(mesh1,mesh2,s1,s2)
3 % mesh=Ujoinmeshes(mesh1,mesh2,side1,side2)
7 % SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator
8 % -------------------------------------------------------------------
9 % Copyright (C) 2004-2006 Carlo de Falco
13 % SECS2D is free software; you can redistribute it and/or modify
14 % it under the terms of the GNU General Public License as published by
15 % the Free Software Foundation; either version 2 of the License, or
16 % (at your option) any later version.
18 % SECS2D is distributed in the hope that it will be useful,
19 % but WITHOUT ANY WARRANTY; without even the implied warranty of
20 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 % GNU General Public License for more details.
23 % You should have received a copy of the GNU General Public License
24 % along with SECS2D; If not, see <http://www.gnu.org/licenses/>.
27 % make sure that the outside world is always
28 % on the same side of the boundary of mesh1
29 [mesh1.e(6:7,:),I] = sort(mesh1.e(6:7,:));
30 for ic=1:size(mesh1.e,2)
31 mesh1.e(1:2,ic) = mesh1.e(I(:,ic),ic);
40 side1 = s1(is);side2 = s2(is);
41 [i,j] = find(mesh1.e(5,:)==side1);
43 [i,j] = find(mesh2.e(5,:)==side2);
44 oldregion(side1) = max(max(mesh2.e(6:7,j)));
48 intnodes1=[mesh1.e(1,j1),mesh1.e(2,j1)];
49 intnodes2=[mesh2.e(1,j2),mesh2.e(2,j2)];
51 intnodes1 = unique(intnodes1);
52 [tmp,I] = sort(mesh1.p(1,intnodes1));
53 intnodes1 = intnodes1(I);
54 [tmp,I] = sort(mesh1.p(2,intnodes1));
55 intnodes1 = intnodes1(I);
57 intnodes2 = unique(intnodes2);
58 [tmp,I] = sort(mesh2.p(1,intnodes2));
59 intnodes2 = intnodes2(I);
60 [tmp,I] = sort(mesh2.p(2,intnodes2));
61 intnodes2 = intnodes2(I);
63 % delete redundant edges
67 indici=[];consecutivi=[];
68 indici = unique(mesh2.e(5,:));
69 consecutivi (indici) = [1:length(indici)]+max(mesh1.e(5,:));
70 mesh2.e(5,:)=consecutivi(mesh2.e(5,:));
73 % change node indices in connectivity matrix
75 indici=[];consecutivi=[];
76 indici = 1:size(mesh2.p,2);
77 offint = setdiff(indici,intnodes2);
78 consecutivi (offint) = [1:length(offint)]+size(mesh1.p,2);
79 consecutivi (intnodes2) = intnodes1;
80 mesh2.e(1:2,:)=consecutivi(mesh2.e(1:2,:));
81 mesh2.t(1:3,:)=consecutivi(mesh2.t(1:3,:));
84 % delete redundant points
85 mesh2.p(:,intnodes2) = [];
88 regions = unique(mesh1.t(4,:));
89 newregions(regions) = 1:length(regions);
90 mesh1.t(4,:) = newregions(mesh1.t(4,:));
93 regions = unique(mesh2.t(4,:));
94 newregions(regions) = [1:length(regions)]+max(mesh1.t(4,:));
95 mesh2.t(4,:) = newregions(mesh2.t(4,:));
96 % set adjacent region numbers in edge structure 2
97 [i,j] = find(mesh2.e(6:7,:));
99 mesh2.e(i,j) = newregions(mesh2.e(i,j));
100 % set adjacent region numbers in edge structure 1
101 mesh1.e(6,j1) = newregions(oldregion(mesh1.e(5,j1)));
103 % make the new p structure
104 mesh.p = [mesh1.p mesh2.p];
105 mesh.e = [mesh1.e mesh2.e];
106 mesh.t = [mesh1.t mesh2.t];
109 % %double check to avoid degenerate triangles
110 % [p,ii,jj]=unique(mesh.p(1:2,:)','rows');
112 % mesh.e(1:2,:)=jj(mesh.e(1:2,:));
113 % mesh.t(1:3,:)=jj(mesh.t(1:3,:));
115 % [ii,jj] = find (mesh.e(1,:)==mesh.e(2,:));
117 % [ii,jj] = find ((mesh.t(1,:)==mesh.t(2,:))|(mesh.t(1,:)==mesh.t(3,:))|(mesh.t(3,:)==mesh.t(2,:)));