octave_packages/secs2d-0.0.8/Utilities/Usubdomains2.m

   1 function [e,t]=Usubdomains2(p,t,rcts,sidelist);
   2
   3 %  [e,t]=Usubdomains(p,t,rcts,sidelist);
   4
   5 e=[];
   6
   7 %%% subdivide domain according to position of
   8 %%% elements' center of mass
   9
  10 % get elements' center of mass coordinates
  11 x = p(1,:);y = p(2,:);
  12 mx = sum(x(t(1:3,:)),1)/3;
  13 my = sum(y(t(1:3,:)),1)/3;
  14
  15 t(4,:) = 1;
  16
  17 % loop over rectangular regions
  18 for ii = 1:size(rcts,1)
  19
  20     % find elements with center of mass in this rectangle
  21     trs = find ((mx>rcts(ii,1))&(mx<rcts(ii,2))&...
  22     (my>rcts(ii,3))&(my<rcts(ii,4)));
  23
  24     % set subdomain number
  25     t(4,trs) = ii+1;
  26
  27 end
  28
  29 % get all element edges
  30 sides = [t([1 2 4],:),t([2 3 4],:),t([3 1 4],:)];
  31 sides(4,:) = 0;
  32 ns = size(sides,2);
  33
  34
  35 % build list of edges using conditions on
  36 % segment vertex coordinates
  37
  38 x1 = p(1,sides(1,:)) ;
  39 x2 = p(1,sides(2,:)) ;
  40 y1 = p(2,sides(1,:)) ;
  41 y2 = p(2,sides(2,:)) ;
  42
  43 e=[];
  44 for icond=1:length(sidelist)
  45
  46 onside  = find( ...
  47 (x1<=sidelist(icond,1)) &...
  48 (x2<=sidelist(icond,2)) &...
  49 (x1>=sidelist(icond,3)) &...
  50 (x2>=sidelist(icond,4)) &...
  51 (y1<=sidelist(icond,5)) &...
  52 (y2<=sidelist(icond,6)) &...
  53 (y1>=sidelist(icond,7)) &...
  54 (y2>=sidelist(icond,8)) );
  55 eonside = unique(sides([1,2],onside)','rows')';
  56 eonside(5,:) = icond;
  57
  58 e=[e,eonside];
  59 end
  60
  61 % set left and right subdomain
  62 for ie = 1:size(e,2)
  63 for it=1:size(t,2)
  64     if(((e(1,ie)==t(1,it))&(e(2,ie)==t(2,it)))|...
  65             ((e(1,ie)==t(2,it))&(e(2,ie)==t(3,it)))|...
  66             ((e(1,ie)==t(3,it))&(e(2,ie)==t(1,it))))
  67             e(6,ie)=t(4,it);
  68     end
  69     if(((e(2,ie)==t(1,it))&(e(1,ie)==t(2,it)))|...
  70             ((e(2,ie)==t(2,it))&(e(1,ie)==t(3,it)))|...
  71             ((e(2,ie)==t(3,it))&(e(1,ie)==t(1,it))))
  72             e(7,ie)=t(4,it);
  73     end
  74 end
  75 end