MATLAB: Using a for loop to change a matrix

circulationequationsfor looppolygontriangles

Hello everyone
Okay, so I have a 40×3 matrix with nodal point IDs. This can then be broken down into x and y coordinates per nodal point. Each line in the 40×3 matrix represents a triangle. I'm trying to make node identification go in a counter-clockwise direction. I've used the equation in the code below to calculate the polygon (triangle) area. For those values that are positive, the direction is counter-clockwise. For negative values it is clockwise. Its roughly 50/50.
What I'm trying to do is use a for loop, so some other method, to make the nodal identification direction change. So instead of having x1 y1, x2 y2, x3 y3, I want x1 y1, x3 y3, x2 y2. I've been mulling over this for a few days, looking online to see if anyone has had similar ideas but I can't find anything. For what I have down below, I keep getting an error. And obviously that wouldn't be the full way of doing it, but it was a way I was going to try.
Any help would be greatly appreciated.
My code:
%%Determing boundary nodes
Tri.triangles=R.ci;
Tri.triX=x(Tri.triangles); Tri.triY=y(Tri.triangles);
x1=Tri.triX(:,1); y1=Tri.triY(:,1);
x2=Tri.triX(:,2); y2=Tri.triY(:,2);
x3=Tri.triX(:,3); y3=Tri.triY(:,3);
% direction of triangle circulation, +ive is counterclockwise
for i=1:length(Tri.triangles);
Tri.area(i)=0.5*((x1(i)*y2(i))-(x2(i)*y3(i))+(x2(i)*y3(i))-(x3(i)*y2(i)));
T_area=Tri.area';
t_neg(i)=T_area(i)<0; T_neg=t_neg';
end
T_neg=double(T_neg);  
bb=x1(T_neg(x1)>0);
Error:
Subscript indices must either be real positive integers or logicals.
T_area=
-0.338828523177654
-0.565042113652453
-0.903958166425582
0.00953685300191864
-0.271408101485577
-0.670282627630513
0.370416248624679
0.260899886256084
-0.225398410111666
-0.00953673134790734
-0.826227724959608
0.826595614606049
-0.418818252452184
0.117110523569863
0.0329695575637743
-0.0329722759779543
-1.57312464073766
-0.658750543370843
0.658923690207303
-0.155193869140930
-0.0884571871720254
0.275831339182332
0.0546458313474432
-0.325210246955976
-0.117114751192275
0.286314578843303
0.494672991335392
0.569307863712311
0.871622048784047
0.121534811332822
1.57387580676004
0.169117909390479
-0.324359477497637
-0.185630145482719
-0.402470884844661
0.402549884282053
-0.286296253907494
0.0546600596280769
0.403549019247294
0.0615939649869688

Best Answer

Given a function that returns the vector of areas given the node and coordinate arrays, then
ix=triArea(nodes,coordinates)<0;          % the rows that fail
nodes(ix,[2:3])=fliplr(nodes(ix,[2:3]));  % swap the offending ids
should do the trick, methinks...
ADDENDUM
OK, am back from meeting and had chance to look into the guts a little more. The area formula as modified in your recent comment still isn't correct; the vector cross formula for area is A=|AB x AC|/2. Algebraically if name the points ABC-->123 this reduces to
AB=[x2−x1,y2−y1]; AC=[x3−x1,y3−y1];  % vectors 
A=[(x2−x1)(y3−y1)−(x3−x1)(y2−y1]/2;  % Area of ΔABC
You can write this in Matlab as a function handle succinctly as:
>> fn_triArea=@(x,y) ((x(:,2)-x(:,1)).*(y(:,3)-y(:,1))-(x(:,3)-x(:,1)).*(y(:,2)-y(:,1)))/2;
In use it's simply
>> A=fn_triArea(x,y)
A =
 1.0e-03 *
  0.3445
  0.2559
  0.2817
  0.2110
  0.3301
  0.2789
>>
For comparison,
>> [A polyarea(x',y').']
ans =
 1.0e-03 *
  0.3445    0.3445
  0.2559    0.2559
  0.2817    0.2817
  0.2110    0.2110
  0.3301    0.3301
  0.2789    0.2789
>>
we note the same values as the builtin polyarea function which gives confirmation we got the formula right. I did a second and third comparison as well of the area of the total rectangle bounding the triangle less the three triangles' outside the inscribed one and the direct cross calculation as well--all were in agreement. I didn't check your algebra but you've still got a transcription or sign error somewhere.
Note the x, y coordinate arrays you've given all result in positive areas, if we pick a couple at random and flip coordinates what happens?
>> ix=randi(length(x),2,1)  % two random points to flip coordinates
ix =
   2
   5
>> x(ix,[2:3])=fliplr(x(ix,[2:3]));  % and flip 'em...
>> y(ix,[2:3])=fliplr(y(ix,[2:3]));
>> [fn_triArea(x,y) A polyarea(x.',y.').']
ans =
 1.0e-03 *
  0.3445    0.3445    0.3445
 -0.2559    0.2559    0.2559
  0.2817    0.2817    0.2817
  0.2110    0.2110    0.2110
 -0.3301    0.3301    0.3301
  0.2789    0.2789    0.2789
>>
Note the two selected triangles now do have negative areas and again compare numerically with the previous results.
Of course, you also need to swap the node numbers to match but if you're storing the full x,y arrays and computing based on them and their row indices instead of actual node ID then the area computation doesn't even need the node; it's just an external bookkeeping addon.
But, there's no indexing error and the logic described will fix your problem implemented correctly.
HTH...
PS: The above results are from your x,y arrays you pasted in in a comment some time earlier that are roughly -6.3 and 56.8 with small differences in values, hence the small areas aren't surprising.
ADDENDUM 2
As one last illustration on the area, the worst possible way to compute anything almost always, but we'll trot it out for further confirmation--
>> det([x(1,:)',y(1,:)',ones(3,1)])/2
ans =
 3.4450e-04
>>
Related Question