% script file tri_select.m
% for a given triangulation and const. vector field a,
% find the 1-inflow and 2-inflow triangles
% last changes: 24.06.03 by Vita
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; close all

a = [1,0.1]; % the given vector field

% load the data files
tri = load('xshape.tri'); % triangles
nod = load('xshape.nod'); % nodes

% number of triangles & nodes
n.tri = length(tri); n.nod = length(nod);

% first, visualize to be sure... :)
figure
trimesh(tri,nod(:,1),nod(:,2),ones(n.nod,1));
view(2); title('the original triangulation')

% compute the inner normals for each triangle
% a new data structures nor1, nor2, nor3
% (for 1-st, 2-nd and 3-rd edge resp.)
% each with 2 columns normal_1 and normal_2

nor1 = zeros(n.tri,2);
nor2 = zeros(n.tri,2);
nor3 = zeros(n.tri,2);

for ctr = 1:n.tri
  % select coordinates of each node
  x1 = nod(tri(ctr,1),1); y1 = nod(tri(ctr,1),2);  
  x2 = nod(tri(ctr,2),1); y2 = nod(tri(ctr,2),2);  
  x3 = nod(tri(ctr,3),1); y3 = nod(tri(ctr,3),2); 
  
  nor1(ctr,1) = y2-y3; nor1(ctr,2) = x3-x2;
  nor2(ctr,1) = y3-y1; nor2(ctr,2) = x1-x3;
  nor3(ctr,1) = y1-y2; nor3(ctr,2) = x2-x1;

  % check if the points are positively oriented
  % for that, check if vectors for edges
  % 1 and 2 are positively oriented: use determinant condition
  % a and b are positively oriented if
  % |a1 b1|
  % |a2 b2| > 0

  D = (x3-x2)*(y1-y3)-(y3-y2)*(x1-x3);  
  if D<0 
    % change the orientation
    disp('directions changed')
    nor1(ctr,:) = -nor1(ctr,:);
    nor2(ctr,:) = -nor2(ctr,:);    
    nor3(ctr,:) = -nor3(ctr,:);  
  end 
end

% compute 1 and 2 inflow edge triangles
% tri_mark = 1 for 1 inflow edge, 2 for 2 inflow edges
tri_mark = zeros(n.tri,1);

for ctr = 1:n.tri
  % compute the scalar products of normals with a
  sp1 = a*nor1(ctr,:)';
  sp2 = a*nor2(ctr,:)';
  sp3 = a*nor3(ctr,:)';
  
  % for 2 inflow edges, two of them are positive,
  % one is negative =>sp1*sp2*sp3<0
  % for 1 inflow edge, two of them are negative,
  % => sp1*sp2*sp3>0
  
  if sp1*sp2*sp3<0
    tri_mark(ctr) = 2;
  else
    tri_mark(ctr) = 1;
  end
end

  
% visualize the results
figure
hold on
for ctr = 1:n.tri
  X = [nod(tri(ctr,1),1);nod(tri(ctr,2),1);nod(tri(ctr,3),1)];
  Y = [nod(tri(ctr,1),2);nod(tri(ctr,2),2);nod(tri(ctr,3),2)];
  fill(X,Y,tri_mark(ctr));
end
colormap([0 1 0;1 0 0]);
colorbar

% plot the vector field
quiver(nod(:,1),nod(:,2),a(1)*ones(n.nod,1), a(2)*ones(n.nod,1));
       
  



  

