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
>>
