# Using a for loop to change a matrix

3 views (last 30 days)
Meghan on 4 Oct 2016
Commented: Meghan on 5 Oct 2016
Hello everyone
Okay, so I have a 40x3 matrix with nodal point IDs. This can then be broken down into x and y coordinates per nodal point. Each line in the 40x3 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

dpb on 4 Oct 2016
Edited: dpb on 4 Oct 2016
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...
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=[x2x1,y2y1]; AC=[x3x1,y3y1]; % vectors
A=[(x2x1)(y3y1)(x3x1)(y2y1]/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.
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
>>
dpb on 5 Oct 2016
Return a list of just the unique nodes (use nodes(:,1) for input) and then either just loop over that list returning the collection for which u(i)==nodes(:,1) and find the max of the code number in each set, returning the optional index variable which will be the index within the set of the desired row. The index array from the unique call will let you relate that row back to the row in the original array.
Or, you can get more sophisticated and take the above outputs from the second unique call and bin them with histc from which you can find the subset of duplicates. accumarray can then be used with an anonymous function including max for each of these groups to do the function of the above "brute force" loop construct. In reality, the loop probably will run faster than the accumarray solution and is simpler to code but by some measures the latter is "more elegant".
Meghan on 5 Oct 2016
Thanks :) I was thinking along the same lines, just wanted to see if you had a simple solution I might have missed. I tend to overthink a lot of things!
I very much appreciate all the help :)