octave_packages/secs2d-0.0.8/Utilities/Ujoinmeshes.m

   1 function mesh=Ujoinmeshes(mesh1,mesh2,s1,s2)
   2
   3 %  mesh=Ujoinmeshes(mesh1,mesh2,side1,side2)
   4
   5 % This file is part of
   6 %
   7 %            SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator
   8 %         -------------------------------------------------------------------
   9 %            Copyright (C) 2004-2006  Carlo de Falco
  10 %
  11 %
  12 %
  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.
  17 %
  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.
  22 %
  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/>.
  25
  26
  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);
  32 end
  33
  34 intnodes1=[];
  35 intnodes2=[];
  36
  37
  38 j1=[];j2=[];
  39 for is=1:length(s1)
  40     side1 = s1(is);side2 = s2(is);
  41     [i,j] = find(mesh1.e(5,:)==side1);
  42     j1=[j1 j];
  43     [i,j] = find(mesh2.e(5,:)==side2);
  44     oldregion(side1) = max(max(mesh2.e(6:7,j)));
  45     j2=[j2 j];
  46 end
  47
  48 intnodes1=[mesh1.e(1,j1),mesh1.e(2,j1)];
  49 intnodes2=[mesh2.e(1,j2),mesh2.e(2,j2)];
  50
  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);
  56
  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);
  62
  63 % delete redundant edges
  64 mesh2.e(:,j2) = [];
  65
  66 % change edge numbers
  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,:));
  71
  72
  73 % change node indices in connectivity matrix
  74 % and edge list
  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,:));
  82
  83
  84 % delete redundant points
  85 mesh2.p(:,intnodes2) = [];
  86
  87 % set region numbers
  88 regions = unique(mesh1.t(4,:));
  89 newregions(regions) = 1:length(regions);
  90 mesh1.t(4,:) = newregions(mesh1.t(4,:));
  91
  92 % set region numbers
  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,:));
  98 i = i+5;
  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)));
 102
 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];
 107
 108 %
 109 % %double check to avoid degenerate triangles
 110 % [p,ii,jj]=unique(mesh.p(1:2,:)','rows');
 111 % mesh.p =p';
 112 % mesh.e(1:2,:)=jj(mesh.e(1:2,:));
 113 % mesh.t(1:3,:)=jj(mesh.t(1:3,:));
 114 %
 115 % [ii,jj] = find (mesh.e(1,:)==mesh.e(2,:));
 116 % mesh.e(:,jj) = [];
 117 % [ii,jj] = find ((mesh.t(1,:)==mesh.t(2,:))|(mesh.t(1,:)==mesh.t(3,:))|(mesh.t(3,:)==mesh.t(2,:)));
 118 % mesh.t(:,jj) = [];
 119
 120