octave_packages/secs2d-0.0.8/Utilities/Umeshproperties.m

   1 function mesh=Umeshproperties(mesh)
   2
   3 %
   4 % mesh=Umeshproperties(mesh)
   5 % precomputes some useful mesh properties
   6 %
   7
   8
   9
  10 % This file is part of
  11 %
  12 %            SECS2D - A 2-D Drift--Diffusion Semiconductor Device Simulator
  13 %         -------------------------------------------------------------------
  14 %            Copyright (C) 2004-2006  Carlo de Falco
  15 %
  16 %
  17 %
  18 %  SECS2D is free software; you can redistribute it and/or modify
  19 %  it under the terms of the GNU General Public License as published by
  20 %  the Free Software Foundation; either version 2 of the License, or
  21 %  (at your option) any later version.
  22 %
  23 %  SECS2D is distributed in the hope that it will be useful,
  24 %  but WITHOUT ANY WARRANTY; without even the implied warranty of
  25 %  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  26 %  GNU General Public License for more details.
  27 %
  28 %  You should have received a copy of the GNU General Public License
  29 %  along with SECS2D; If not, see <http://www.gnu.org/licenses/>.
  30
  31
  32
  33
  34 weight = [1/3 1/3 1/3];
  35 areakk = 1/2;
  36 Nelements = Ucolumns(mesh.t);
  37
  38 jac([1,2],:) = [mesh.p(1,mesh.t(2,:))-mesh.p(1,mesh.t(1,:));
  39             mesh.p(1,mesh.t(3,:))-mesh.p(1,mesh.t(1,:))];
  40 jac([3,4],:) = [mesh.p(2,mesh.t(2,:))-mesh.p(2,mesh.t(1,:));
  41             mesh.p(2,mesh.t(3,:))-mesh.p(2,mesh.t(1,:))];
  42 jacdet = jac(1,:).*jac(4,:)-jac(2,:).*jac(3,:);
  43
  44 degen=find(jacdet <= 0);
  45 if ~isempty(degen)
  46   fprintf(1,'invalid mesh element:  %d  fixing...\n',degen);
  47   mesh.t(1:3,degen) = mesh.t([2,1,3],degen);
  48   jac([1,2],degen) = [mesh.p(1,mesh.t(2,degen))-mesh.p(1,mesh.t(1,degen));
  49                       mesh.p(1,mesh.t(3,degen))-mesh.p(1,mesh.t(1,degen))];
  50   jac([3,4],degen) = [mesh.p(2,mesh.t(2,degen))-mesh.p(2,mesh.t(1,degen));
  51                       mesh.p(2,mesh.t(3,degen))-mesh.p(2,mesh.t(1,degen))];
  52   jacdet(degen) = jac(1,degen).*jac(4,degen)-jac(2,degen).*jac(3,degen);
  53 end
  54
  55 for inode = 1:3
  56   mesh.wjacdet(inode,:) = areakk .* jacdet .* weight(inode);
  57 end
  58
  59 mesh.shp     = eye(3);
  60
  61 x0 = mesh.p(1,mesh.t(1,:));
  62 y0 = mesh.p(2,mesh.t(1,:));
  63 x1 = mesh.p(1,mesh.t(2,:));
  64 y1 = mesh.p(2,mesh.t(2,:));
  65 x2 = mesh.p(1,mesh.t(3,:));
  66 y2 = mesh.p(2,mesh.t(3,:));
  67
  68 denom = (-(x1.*y0) + x2.*y0 + x0.*y1 - x2.*y1 - x0.*y2 + x1.*y2);
  69 mesh.shg(1,1,:)  =  (y1 - y2)./denom;
  70 mesh.shg(2,1,:)  = -(x1 - x2)./denom;
  71 mesh.shg(1,2,:)  = -(y0 - y2)./denom;
  72 mesh.shg(2,2,:)  =  (x0 - x2)./denom;
  73 mesh.shg(1,3,:)  =  (y0 - y1)./denom;
  74 mesh.shg(2,3,:)  = -(x0 - x1)./denom;
  75
  76