1 function [e,t]=Usubdomains2(p,t,rcts,sidelist);
3 % [e,t]=Usubdomains(p,t,rcts,sidelist);
7 %%% subdivide domain according to position of
8 %%% elements' center of mass
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;
17 % loop over rectangular regions
18 for ii = 1:size(rcts,1)
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)));
24 % set subdomain number
29 % get all element edges
30 sides = [t([1 2 4],:),t([2 3 4],:),t([3 1 4],:)];
35 % build list of edges using conditions on
36 % segment vertex coordinates
38 x1 = p(1,sides(1,:)) ;
39 x2 = p(1,sides(2,:)) ;
40 y1 = p(2,sides(1,:)) ;
41 y2 = p(2,sides(2,:)) ;
44 for icond=1:length(sidelist)
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')';
61 % set left and right subdomain
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))))
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))))